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constitutes  about  2%  of  the  incident  ambient  flux;  this  ratio  will  be  nearly  constant  for  LEO 
altitudes.  The  value  of  this  flux  ratio  is  shown  to  be  dependent  on  the  molecular  collision  model;  it 
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1.        INTRODUCTION 

This  presentation  is  part  of  a  study  on  the  gas  dynamics  of  ring-symmetric  exhaust  plumes  in 
space,  conducted  at  the  Naval  Postgraduate  School  in  Monterey.  A  ring- symmetric  jet  has  zero 
thrust,  which  makes  it  suitable  as  an  exhaust  configuration  for  various  open  loop  power  plants 
designed  to  produce  high  power  for  relatively  short  durations.  One  such  system  is  an  envisioned 
space-based  chemical  laser,  shown  schematically  in  Fig.  1-1.  In  the  case  of  a  chemical  laser,  a  ring- 
symmetric  configuration  would  also  enable  the  laser  radiation  to  emerge  in  the  form  of  an 
axisymmetric  beam. 

The  exhaust  nozzle  should  be  designed  to  bring  the  outgoing  flow  to  a  supersonic  speed  at  the 
nozzle  exit  surface.  The  near  field  of  a  free  jet  is  then  composed  of  an  inner  core  bounded  by  a  pair 
of  ring- symmetric  rarefaction  fans  centered  at  the  nozzle  lips  (Fig.  1-1).  Beyond  the  limiting 
characteristic  surface  of  the  centered  rarefaction  waves  (CRW),  a  near-vacuum  condition  prevails. 
For  the  purpose  of  continuum  gas  dynamic  analysis,  we  assume  it  is  a  perfect  vacuum. 

An  earth  orbiting  vehicle  is  subject  to  an  oncoming  stream  of  ambient  molecules  at  a  speed  of 
UA*8  (km/sec),  in  a  direction  depending  upon  its  orientation  relative  to  the  orbital  velocity  vector. 
This  speed  is  sufficiently  high  to  cause  backscattering  of  exhaust  molecules  (see  schematic  description 
in  Fig.  1-2)  moving  at  speeds  appropriate  to  chemical  combustion  (about  2  to  4  km/s).  However, 
large  exhaust  plumes,  having  achieved  stationary  flow,  may  be  sufficiently  dense  at  their  outer  fringes 
to  effectively  trap  and  entrain  all  oncoming  ambient  molecules.  Thus,  ambient  scattering  may  be 
significant  only  in  selected  ranges  of  attitude  angles,  at  which  ambient  molecules  can  reach  the 
vicinity  of  the  spacecraft  by  traveling  almost  collisionlessly  through  cavitation  regions.  Exhaust 
molecules  that  may  be  "candidates"  for  ambient  scattering  will  hence  come  from  plume  segments 
flanked  by  cavitation  regions.  The  contribution  of  ambient  scattering  to  contamination  will  thus  be 
highly  dependent  upon  spacecraft  geometry  and  orientation.  This  may  well  affect  spacecraft  design 
and  operating  procedures. 

The  purpose  of  this  report  is  to  present  a  first-collision  model  for  estimating  the  flux  of  exhaust 
molecules  backscattered  from  the  fringes  of  the  plume  by  ambient  molecules,  along  with  results  of 
sample  flux  computations  performed  on  a  typical  HF/DF  laser  exhaust  configuration.  The  flow  field 
throughout  the  plume  is  assumed  to  be  governed  by  the  equations  of  continuum  gas  dynamics.  In 
principle,  the  How  could  be  obtained  by  solving  the  governing  equations,  i.e..  the  equations  for 
stationary  isentropic  flow  in  two-dimensional  axisymmetric  coordinates.    In  practice,  this  is  normally 


accomplished  by  integrating  the  flow  equations  in  characteristic  form,  using  some  finite  difference 
scheme  (method-of-characteristics).  We  have  performed  such  computations,  but  given  the  complexity 
of  applying  them  to  the  subsequent  integration  of  ambient  scattering  flux  (due  to  the  need  for  two- 
dimensional  interpolations  from  an  irregular  solution  grid),  we  opted  for  a  different  alternative  :  a 
closed-form  approximation  to  the  ring- symmetric  CRW,  based  on  an  analytic  expression  for  flow 
variables  along  characteristic  lines  that  fan  out  from  the  nozzle  lip. 

The  plan  of  this  report  is  as  follows.  In  Ch.  2  we  outline  the  approximation  to  the  ring- symmetric 
CRW  and  present  some  computation  results  that  demonstrate  its  accuracy.  In  Ch.  3  we  describe  the 
first-collision  model  and  the  3-D  spatial  integration  scheme  for  computing  the  flux  arriving  at  the 
cylindrical  spacecraft.  In  Ch.  4  some  results  of  backscattered  flux  of  corrosive  molecules  (HF  +  DF), 
showing  flux  variation  with  target  point  location  (Xs)  and  attitude  angles  (\yA,  <pA)  are  presented.  In 
Ch.  5  we  take  up  the  subject  of  spacecraft  charging,  using  results  of  ambient  scattering  to  assess  the 
effect  of  laser  exhaust  on  spacecraft  charging.  This  is  followed  by  concluding  remarks  in  Ch.  6  and  a 
list  of  references  in  Ch.  7.  A  concise  description  of  the  flux  computation  code  "AMB"  is  given  in 
Appendix  A,  followed  by  the  code  listing. 
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Figure    1-1.      Ring-Symmetric  HF  DF  Laser  Exhaust  Plume. 
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Figure    1-2.      Schematic  Description  of  Ambient  Scattering.    The  Cavitation  Region  is  Bounded  by 
Lines  CA  and  CB. 


2.        COMPUTATION   OF  THE   PLUME   FLOW  FIELD 

Most  ambient  molecules  entering  the  CRW  that  flanks  the  exhaust  plume  are  stopped  within 
several  mean  free  paths  from  their  point  of  entry.  A  quantitative  estimate  of  ambient  back-scattering 
would  thus  depend  on  the  flow  field  at  the  outer  (hypersonic)  fringes  of  the  lip-centered  CRW.  Even 
though  the  flow  in  those  regions  is  generally  past  the  surface  of  continuum  breakdown,  the  density 
there  is  reasonably  well  approximated  by  the  continuum  flow  field,  as  demonstrated  by  Bird's  Monte- 
Carlo  simulation  of  a  Prandtl-Meyer  expansion  to  vacuum  [1]  .  The  evaluation  of  ambient  scattering 
thus  calls  for  an  ancillary  computational  procedure  capable  of  rendering  the  continuum  flow  field  at  a 
large  number  of  points  in  the  ring- symmetric  CRW  of  an  exhaust  plume.  This  method  was  described 
in  a  recent  report  [2]  .  Here  we  just  outline  the  key  ideas  and  main  results  of  this  approximation 
method. 

Our  analytic  approximation  to  a  ring- symmetric  CRW  is  formulated  as  follows.  In  a  planar  CRW 
(Prandtl-Meyer  flow)  all  flow  variables  are  uniform  along  the  characteristic  lines  that  fan  out  from  the 
corner  (we  assume  they  are  the  C  family).  In  the  ring- symmetric  case  the  flow  near  the  corner 
approaches  asymptotically  a  corresponding  planar  CRW  flow,  which  we  term  the  associate  CRW. 
However,  the  gradients  along  C  characteristics  at  the  corner  of  a  ring- symmetric  CRW  do  not 
vanish  as  in  a  planar  CRW.  The  key  idea  is  thus  :  evaluate  flow  gradients  in  C  directions  at  the 
corner,  then  use  them  to  extrapolate  the  associate  CRW  along  C  lines  to  a  finite  distance  from 
the  corner.  The  extrapolation  is  a  nonlinear  function  of  the  radial  coordinate  y  ,  chosen  so  that  the 
ensuing  expression  conforms  exactly  to  the  flow  at  the  leading  (exit)  characteristic  C  (Pj)  . 
Omitting  all  details  of  the  analysis,  the  resulting  approximation  is  presented  as  the  following  power- 
law  : 

f(a.P)  =  f(0,P)  [y(a,P)/y(0,P)]  (2-1) 

where  f  is  the  streamtube  area  ratio  for  isentropic  flows  (  f  =  1  at  a  sonic  point),  P  is  the  Mach 
number  of  a  particular  characteristic  line  at  the  corner,  a  is  a  coordinate  along  the  C  (p) 
characteristic  line  (  a  =  0  at  the  corner),  and  y  is  the  radial  coordinate  of  a  point  on  the 
characteristic  line  C  (P)  .  The  Mach  number  at  point  (a,p)  is  readily  determined  from  f(a.p)  using 
the  standard  relation  between  area  ratio  and  Mach  number  [3]  .  A  closed-form  expression  for  6(0. P) 
was  developed  but  is  not  given  here;  instead,  this  function  is  shown  in  Fig.  2-1.  We  note  that  6 
approaches  the  asymptotic  value  of  2/(3-y)  as  P  increases  to  infinity,  and  that  generally 
1  <  5(0, P)  <  2  so  that  streamtubes  diverge  at  a  rate  intermediate  between  that  of  cylindrical  and 
spherical  expansion  flows. 


Clearly,  in  an  isentropic  flow  all  thermodynamic  variables,  and  in  particular  density,  can  be 
evaluated  from  f .  This  approximation  is  readily  applied  to  the  hypersonic  portions  of  a  ring- 
symmetric  CRW  since  it  turns  out  that  characteristic  lines  are  nearly  straight  there,  which  means  that 
the  characteristic  line  C  (P)  passing  through  a  given  point  can  be  readily  determined.  As  a 
demonstration  of  the  degree  of  accuracy  obtainable  from  this  approximation,  we  show  in  Fig.  2-2  the 
variation  of  Mach  number  along  a  characteristic  line  in  the  ring- symmetric  CRW,  compared  with  an 
accurate  method-of-characteristics  computation.  This  comparison  demonstrates  that  the  analytic 
approximation  is  reasonably  accurate  to  nearly  ten  corner-radii  away  from  the  corner. 
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Figure  2-1.      Power  6(0,0)  for  the  power-law  Approximation. 
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Figure  2-2.      Variation  of  Mach  Number  along  Characteristic  Line  p=  13. 


3.        AMBIENT  SCATTERING 

When  a  rocket  or  laser  exhaust  is  released  into  space  from  an  earth-orbiting  spacecraft,  it 
encounters  an  oncoming  stream  of  ambient  molecules  flowing  at  the  orbital  speed  of  UA*8  (km/sec). 
At  altitudes  higher  than  200  (km),  the  air/air  mean  free  path  exceeds  250  (m),  so  that  it  is 
considerably  larger  than  almost  any  spacecraft.  Consequently,  ambient  molecules  would  hardly  be 
subjected  to  a  shock  transition  prior  to  their  impact  at  the  spacecraft  or  exhaust  plume.  In  this 
chapter  we  describe  the  formulation  of  the  first-collision  model  in  Section  3.1  and  then  proceed  to 
present  the  derivation  of  the  flux  integration  scheme  for  hard-sphere  collisions  in  Section  3.2. 


3.1      First  Collision  Model 

The  highest  ambient  number  density  that  we  consider  for  earth- orbiting  spacecrafts  is 
NA=  1  x  1016  (m"3),  which  roughly  corresponds  to  Sunspot  Maximum  at  200  (km)  [4]  .  The  typical 
laser  exhaust  (Table  4-1)  would  reach  a  number  density  of  about  2  x  1019  (m"3)  at  the  very  high  Mach 
number  of  30.  Hence,  ambient  flux  constitutes  just  a  slight  perturbation  to  the  near-field  portion  of 
a  typical  laser  exhaust  plume.  Obviously,  ambient  molecules  that  penetrate  the  plume,  would 
subsequently  be  entrained  by  the  main  flow.  But  how  far  do  they  penetrate?  And  would  exhaust 
molecules  scattered  by  them  reach  the  spacecraft?  In  seeking  answers  to  these  questions,  we  are  led 
to  some  interesting  observations  concerning  ambient  scattering. 

Consider  the  HF  laser  depicted  in  Fig.  1-1.  The  spacecraft  diameter  is  5  (m)  and  the  centrally 
located  ring- symmetric  nozzle  is  1  (m)  wide.  Typical  operating  conditions  (Table  4-1)  are  assumed. 
They  are  based  on  some  experimental  HF/DF  laser  studies  conducted  at  TRW  [5.6]  .  Suppose  that 
the  spacecraft  axis  is  normal  to  the  orbital  velocity  vector  (normal  incidence).  Let  the  plane  of 
incidence  be  the  plane  defined  by  the  intersection  of  the  spacecraft  axis  with  the  orbital  velocity 
vector.  The  probability  that  an  ambient  molecule  traveling  in  the  plane  of  incidence  would  reach  the 
spacecraft  collisionlessly  is  exp(  —  nj,  where  rj  is  its  expected  number  of  collisions  with  exhaust 
molecules.  We  define  the  number  y\  as  "molecular  thickness",  in  analogy  to  "optical  thickness".  So  in 
order  to  determine  the  extent  to  which  ambient  molecules  at  normal  incidence  reach  the  spacecraft, 
we  seek  the  distribution  of  radial  molecular  thickness  as  function  of  distance  from  the  spacecraft  mid- 
plane  (normal  to  axis  at  its  midpoint). 


For  this  purpose  we  computed  the  ring-symmetric  exhaust  flow  field,  using  a  semi-inverse 
marching  characteristics  scheme  [7]  .  The  marching  was  in  the  radial  direction,  starting  with  uniform 
flow  at  the  nozzle  exit;  the  computation  was  carried  on  until  it  became  evident  that  even  at  a 
distance  of  20  (m)  from  the  mid-plane,  the  radial  molecular  thickness  was  well  over  40.  The  entire 
spacecraft  was  thus  shielded  from  any  ambient  scattering  at  (or  near)  normal  incidence.  This 
shielding  effect  has  two  significant  implications  which  we  discuss  briefly  below. 

(a)  It  is  present  only  during  stationary  exhaust  flow.  At  startup  and  shutdown  phases,  ambient 
scattering  may  be  substantial  even  at  normal  incidence. 

(b)  During  the  stationary  phase,  ambient  scattering  is  substantial  only  at  attitude  angles  that  enable 
ambient  molecules  to  reach  the  vicinity  of  the  plume  by  traveling  through  "molecularly  thin'' 
cavitation  regions  that  flank  the  plume.  We  thus  anticipate  a  decisive  dependence  of  ambient 
scattering  on  attitude  variations,  whenever  those  variations  steer  the  spacecraft  into  or  out  of  a 
shielded  posture. 

As  a  first  attempt  at  a  quantitative  estimate  of  ambient  scattering  flux,  we  have  formulated  a 
simple  first-collision  model  of  this  effect.  In  the  sequel  we  present  an  outline  of  the  model,  along  with 
some  sample  results  evaluated  for  an  HF  laser  configuration  identical  to  that  considered  for  the 
shielding  effect  mentioned  above. 

The  basic  idea  is  the  following.  Ambient  molecules  entering  an  exhaust  plume,  require  several 
collisions  to  become  fully  "accommodated"  with  the  main  flow  (i.e.,  to  be  entrained  by  the  main  How 
at  the  prevailing  flow  velocity  and  temperature).  One  may  reasonably  approximate  this  process  by 
considering  just  one  collision  -  the  first. 

With  the  help  of  some  additional  assumptions,  we  were  able  to  derive  a  closed  form  expression  for 
the  flux  of  exhaust  molecules  that  arrive  at  the  spacecraft  following  a  first  collision  with  an  ambient 
molecule.   The  main  assumptions  of  this  model  are  : 

(1)  FIRST  COLLISIONS:  Only  first  collisions  for  either  ambient  or  exhaust  molecules  are 
considered.  Hard-spheres  elastic  collisions  are  assumed.  Upon  a  second  collision  of  either  an 
ambient  or  an  exhaust  molecule,  it  is  considered  "lost"  (i.e.,  it  joins  the  main  flow).  Collisions 
of  ambient  molecules  with  spacecraft  surfaces  are  ignored.  Ambient  molecules  are  assumed  to 
traverse  cavitation  regions  collisionlessly. 


(2)  COLD  FLOW:  The  oncoming  ambient  air  flow  is  deemed  "cold";  i.e.,  all  molecules  move  at 
the  uniform  orbital  velocity.  The  same  "cold"  assumption  is  applied  to  the  exhaust  flow,  since 
most  ambient  scattering  takes  place  at  plume  regions  of  very  high  Mach  numbers  (well  over  10, 
in  the  present  case). 

(3)  CRW  Flow  Field:  ring- symmetric  CRW  flow  field  is  determined  from  the  power-law 
approximation  described  in  Ch.  2  above.  This  approximation  approaches  Prandtl- Meyer  flow 
at  points  whose  distance  from  the  nozzle  lip  is  much  smaller  than  the  spacecraft  radius. 

Based  on  these  assumptions,  ambient  scattering  is  represented  as  a  source  term  for  side-scattered 
exhaust  molecules,  distributed  throughout  the  lip-centered  rarefaction  fan.  The  total  flux  arriving  at  a 
specified  point  on  the  cylindrical  spacecraft  is  readily  computed  by  integrating  numerically  that  source 
distribution  over  the  entire  ring-fan. 

The  highlights  of  the  spatial  integration  scheme  (Fig.  3-1)  are  as  follows.  The  limiting 
characteristic  surface  (M=  oo)  of  the  ring- symmetric  CRW  is  divided  into  surface  elements  formed  by 
dividing  the  surface  into  a  set  of  ring-strips  which  are  subdivided  in  the  circumferential  (azimuthal) 
direction  (q>)  into  surface  elements.  The  line-of-sight  {ft)  from  the  "target  point"  on  the  spacecraft  to 
the  center  of  each  surface  element  is  extended  into  the  ring- symmetric  CRW,  and  flux  integration 
using  the  first-collision  source  term  with  appropriate  weight  factors  is  performed  along  this  line  until 
convergence  is  attained.  Contributions  from  each  surface  element  are  summed,  taking  care  to 
disregard  portions  of  the  ring-symmetric  CRW  that  are  shadowed  by  the  cylindrical  spacecraft  (either 
the  line-of-sight  or  the  trajectory  of  oncoming  ambient  molecules  may  be  shadowed).  Some  further 
details  of  the  flux  integration  scheme  and  hard-spheres  collisions  are  provided  in  Section  3.2  below. 


3.2      Flux  Integration  Scheme 

The  description  of  the  first  collision  model  is  hereby  supplemented  with  an  outline  of  the 
expressions  used  in  the  flux  integration  and  their  derivation.  The  integration  scheme  for  flux  arriving 
at  point  X§  on  the  spacecraft  is  depicted  in  Fig.  3-1.  Note  that  only  the  plane  of  incidence  is  shown 
in  Fig.  3-1;  at  other  azimuth  angles  the  geometry  is  not  co-planar,  so  3-D  geometrical  expressions  are 
used  to  get  the  coordinates  (\j/.(p  and  radial  distancete  (y2  +  z2)^2)  from  £1  and  S;  the  derivation  of 
these  geometrical  relations  is  straightforward,  so  that  we  omit  these  details  in  the  present  report.  The 
total  number  flux  Qj(Xs)  of  i  exhaust  molecules  arriving  at  point  Xs  is  given  by  the  following 
expression  : 
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Qj(Xs)  =  J  d312  cosas  I  J  dS  (Tik  h.  N(S)  hR  NA  |  U(S)  -  UA  |  expl  -  nk(S)]  P.  (S,  -  ft)  exp[  -  n-  (S)] 

K  o 

t(S) 

Tik(S)  =  2  J  df  (Tik  h.  N(t')  I  U(f)  -  UA I  / 1 UA I  (3-1) 

1  o 

s  ^ 

nik(S)  =  s  J  dS'  (Tjj  h.  N(so  |  uik(S)  -  u(so  I  / 1  uik(S) 


( )j    ( )•    -    Exhaust  species  ( )k    -    Ambient  species 


These  expressions  are  interpreted  as  follows.  The  collision  depicted  in  Fig.  3-1  is  between  exhaust 
molecule  rrij  and  ambient  molecule  mk.  The  exhaust  molar  fractions  h-  and  ambient  molar  fractions 
hk  are  assumed  uniformly  constant,  and  so  are  the  ambient  velocity  UA  and  number  density  N,.  The 
exhaust  velocity  U(S)  and  number  density  N(S)  are  function  of  the  location  in  the  flow  field  defined 
by  12  and  S.  These  flow  variables  are  computed  by  first  evaluating  the  coordinates  of  point  H.S 
(Fig.  3-1)  in  the  ring-symmetric  CRW  from  the  3-D  geometry,  and  then  employing  the  power-law 
approximation  outlined  in  Ch.  2  above,  to  get  all  flow  variables  for  a  ring- symmetric  CRW.  In  this 
computation  we  exploit  the  fact  that  characteristic  lines  fanning  out  from  the  nozzle  lip  are  nearly 
straight  lines  at  the  low  pressure  side  of  the  ring-symmetric  CRW. 

The  Q.  integration  is  performed  numerically  according  to  the  scheme  outlined  in  Section  3.1  above, 
as  a  summation  over  elements  of  solid  angle  (Afl)  subtended  by  area  elements  on  the  limiting 
characteristic  cone  (\|/  =  \j/f). 

The  S  integration  is  considerably  more  complex.  The  integrand  for  this  integration  is  derived  as 
follows.  Denote  by  L  the  line-of-sight  distance  between  point  X  and  fan  point  n,S-  A  volume 
element  at  the  fan  point  is  given  by  Av  =  L2  AS  A3H.  The  number  of  ik  pair  collisions  in  Av  per  unit 
time  is  <Tjk  h.  N(S)  hk  N^  |U(S)-  uj  exp(-  \{S)\  Av  ,  where  \(S)  denotes  the  expected  number  of 
collisions  of  ambient  molecule  k.  with  any  exhaust  molecule,  between  its  point  of  entry  into  the  plume 
and  point  H.S.  We  now  multiply  this  term  by  exp[  -  rjik(S)I  which  is  the  probability  that  exhaust 
molecule  i  scattered  by  ambient  molecule  k  would  travel  from  point  H.S  to  point  X.  collisionlessly, 
where  i\ik{S)  is  the  expected  number  of  collisions  for  this  path  segment.  (Note  that  in  Eq.  (3-1)  the 
summation  in  the  expression  for  Tlik(S)  is  over  all  exhaust  species  j). 
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The  final  step  in  constructing  the  integrand  for  the  S  integration  involves  the  post-collision 
directional  distribution  function  P;k(S,  —  12),  whose  derivation  will  be  given  in  the  sequel.  We 
multiply  the  integrand  by  P,k(S,  — 12)  A  12c  which  is  the  fraction  of  i  exhaust  molecules  scattered  by  k. 
ambient  molecules  into  a  solid  angle  element  A312c  about  the  unit  vector  — 12.  Considering  the  flux 
arriving  at  a  surface  area  element  AAs  around  point  X$,  the  solid  angle  element  subtended  by  AAs  is 
A312  =  AAs  cosas  /L2.  Eq.  (3-1)  for  Q[(Xs)  now  follows  upon  dividing  the  resulting  expression  by 
AA  .  thus  referring  the  arriving  flux  to  a  unit  area  at  the  point  of  arrival  X  . 

Numerically,  the  S  integration  was  performed  using  the  classical  Runge-Kutta  scheme  (fourth 
order).  The  integration  for  n,k(S)  and  \(S)  has  to  be  repeated  at  each  point  S.  We  found 
reasonable  convergence  with  4  points  in  the  \(S)  integration  and  6  points  in  the  azimuth  integration. 
The  S  integration  was  terminated  when  convergence  was  attained  (this  is  the  meaning  of  the  upper 
limit  oo  in  the  S  integral  in  Eq.  (3-1)).  The  summation  over  new  strips  on  the  limiting  cone  (V|/  =  vj/f) 
was  also  terminated  upon  convergence.  The  CPU  time  consumed  per  target  point  was  about  100 
(sec)  on  IBM  3033  mainframe. 

We  now  take  up  the  derivation  of  an  expression  for  the  post-collision  directional  distribution 
function  Pik(S,  —  12),  which  we  denote  hereafter  as  P(  — 12).  We  adopt  the  pair-collision  notation 
presented  in  Fig.  3-2  for  the  hard-sphere  collision  analysis. 

As  a  consequence  of  conservation  of  momentum  and  energy  (elastic  collisions),  the  center-of-mass 
velocitv  C  and  the  magnitude  of  the  relative  velocitv  C  are  unchanged  bv  the  collision  [81  .  The 
post-collision  velocities  are  given  by  : 

C*  =  C    +  n  c»  c*  =  cm  -  n,C* 

I  m         ~i    r  2  m  ~1     r 

r   —       1  2  r  ~~       1  2 

(3-2) 
Hj  =  mll(ml  +  m2)  H2  =  m2/(m1+m2) 

cB-*ci  +  ^c2  |c;|  =  |cr| 

The  only  free  parameter  in  the  expressions  for  post-collision  velocities  is  the  orientation  of  the 
post-collision  relative  velocity  C  .  This  orientation  is  uniformly  likely  to  be  in  any  direction  in  space 
when  hard-spheres  collision  is  assumed  [8]  ,  as  represented  by  the  spherical  scattering  envelope  in 
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Fig.  3-3.    The  probability  of  obtaining  C*  in  solid  angle  element  AJH  about  —fl  (Fig.  3-3)  is  given 
by: 

P(-3)      =  (l/47r|n2Cr|2)  (aa/a33)  =  (l/47t|cos6|)  (|cf|2/|n2Cr|2)  (3-3) 

where  A  A  is  an  area  element  on  the  scattering  envelope,  whose  projection  on  a  plane  normal  to  Q.  is 
coso  I .  We  note  that  the  origin  of  Cm  in  Fig.  3-3  is  external  to  the  scattering  envelope,  resulting 
in  two  possible  scattering  elements  on  the  sphere.  In  all  the  cases  that  we  computed,  however  (see 
Ch.  4  below),  that  point  was  found  to  be  always  internal,  so  that  there  was  only  a  single  scattering 
solution  with  post-collision  velocity  Ujk(S)  pointing  at  the  spacecraft  for  any  ik  pair  collision. 
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Figure   3-1.      Incidence-Plane  Description  of  Flux  Integration  Scheme. 
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Figure   3-2.      Hard-Spheres  Collision  Notation. 
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Figure   3-3.      Scattering  Envelope  for  Hard-Spheres  Collision. 
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4.        RESULTS   AND   DISCUSSION 

We  performed  several  computations  of  return  flux  generated  by  ambient  scattering,  aimed  at 
demonstrating  the  expected  flux  level  and  its  variation  with  spacecraft  target  point  and  orbital 
attitude  angles.  In  all  these  computations  we  assumed  that  the  exhaust  flow  is  as  in  the  typical 
HF/DF  laser  case  (Table  4-1  below),  and  that  the  ambient  density  and  velocity  are  NA=  1  x  1016 
(molecules/m3)  and  UA  =  8  (km/sec).  As  an  approximation  we  further  assumed  that  the  sole  ambient 
species  is  molecular  nitrogen  (molecular  weight  WA  =  28)  and  that  all  binary  collision  cross-sections 
are  uniformly  given  by  <J=7tD  ,  where  D  is  the  molecular  diameter  (Table  4-1).  In  each  computation 
we  evaluated  the  combined  HF  +  DF  flux  by  assuming  that  the  molar  fraction  of  DF  is  zero  and  the 
molar  fraction  of  HF  is  the  combined  value  for  both  species  (Table  4-1)  :  .091  +  .135  =  .226  .  This  is 
justified  by  the  relatively  small  difference  in  molecular  weight  (just  5%)  between  these  two  species. 

Three  sets  of  flux  computation  were  performed  as  follows : 

(a)  Incidence-plane  (cpA  =  0)  target  points  at  various  distances  from  the  nozzle  lip  (Xs  =  .l  to 
X$=  10  (m)),  and  at  constant  incidence  angle  (iyA  =  20°).  The  results  are  shown  in  Fig.  4-1.  We 
observe  that  the  flux  is  fairly  insensitive  to  Xs.  Also  shown  in  Fig.  4-1  are  flux  computations 
where  the  ring- symmetric  CRW  flow  is  approximated  as  a  planar  CRW  (Prandtl- Meyer  How), 
rather  than  the  power-law  as  in  Eq.  (2-1)  above.  The  planar  case  exhibits  a  somewhat  higher 
flux,  particularly  at  large  X.. 

(b)  Incidence-plane  (cp^=0)  target  points  at  Xs=  1  (m)  and  at  various  incidence  angles  (vy^=0  to 
i}/A=40o).  A  polar  representation  of  the  results  is  given  in  Fig.  4-2.  Note  the  sharp  decrease  in 
flux  as  the  incidence  angle  \|/A  approaches  the  plume  limiting  angle  \|/f=41°. 

(c)  Azimuth  angle  variation  (cp  v=  0  to  (pA=  180°)  at  a  constant  location  (Xs=  1  (m))  and  at  a 

constant  angle  of  incidence  (\j/A  =  20°).  A  polar  representation  of  the  results  is  given  in  Fig.  4-3. 
Observe  that  flux  becomes  sensitive  to  azimuth  angle  (p^  only  past  <p^  =  90°,  where  shadowing 
by  the  cylindrical  spacecraft  becomes  increasingly  dominant. 

In  addition  to  return  flux  we  also  computed  the  rms  velocity  of  the  arriving  molecules.  For  the 
target  points  in  group  (a),  the  rms  velocity  varied  between  6000  and  6600  (m/sec)  (the  higher  velocity 
at  smaller  X  ),  which  corresponds  to  a  kinetic  energy  of  about  4  (ev)  per  molecule  (HF). 
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The  maximum  return  flux  arriving  at  the  spacecraft  is  about  0.15  x  1019  (molecules/m2sec),  which 
corresponds  to  a  surface  deposition  rate  of  about  300  monolayers  (HF  +  DF)  per  hour.  This  level  of 
contaminating  flux  may  seem  to  be  not  outright  negligible;  however,  since  return  flux  is  proportional 
to  ambient  density,  it  will  be  scaled  down  considerably  at  higher  altitudes  (and  lower  ambient 
densities). 

We  observe  that  the  maximum  return  flux  constitutes  a  fraction  of  about  2%  of  the  incident 
ambient  flux.  This  return  flux  ratio  is  roughly  maintained  at  almost  all  target  points  and  attitude 
angles  in  groups  (a),  (b)  and  (c).  The  only  exceptions  are  incidence  angles  near  the  limiting  cone 
(y  =  \j/f)  or  at  azimuth  angles  (pA  >  125°  where  shadowing  becomes  dominant.  This  observation  is 
interpreted  as  follows. 

Consider  the  total  solid  angle  subtended  by  the  limiting  cone  (considered  to  be  infinitely  extended 
in  the  axial  direction)  as  viewed  from  a  target  point  (for  all  lines-of-sight  H  pointing  outward  of  the 
cylindrical  spacecraft  surface).  It  is  independent  of  target  location  due  to  the  "self- similar"  geometry. 
During  each  flux  computation,  we  also  evaluated  the  total  solid  angle  subtended  by  that  segment  of 
the  cone  over  which  the  flux  integration  was  actually  performed  (see  Section  3.2).  It  was  found  out 
that  for  all  but  the  "shadowed"  cases  (cpx  >  125°),  this  solid  angle  constituted  a  fraction  of  86  ±  1% 
of  the  solid  angle  subtended  by  the  infinite  cone.  We  interpret  this  result  as  a  hint  that  geometrical 
"view  factors"  arising  in  the  course  of  the  flux  integration,  are  not  the  dominant  factor  in  determining 
the  2%  level  of  flux  ratio.   What  then  are  the  dominant  factors? 

For  a  possible  explanation  we  turn  to  the  flux  integration  scheme  presented  in  Section  3.2.  The 
flux  ratio  is  obtained  upon  dividing  the  integrand  in  Eq.  (3-1)  by  NAUA  and  setting  hR=  1  (since  we 
assume  a  single  species  air).  The  major  factors  in  the  flux  ratio  integrand  appear  to  be  the  no- 
collision  probabilities  exp(  —  Hj^S)]  and  exp[  -  Hk(S)|,  and  the  post-collision  directional  distribution 
function  P;k(S,  — H).  The  flux-averaged  values  of  these  functions  in  the  group  (a)  computations  were 
found  to  be  as  follows  :  Pik(S,-n)  =  .09  to  .10  ,  r|ik(S)  =  .42  to  .54  and  \(S)=  .35  to  .47.  The 
flux-averaged  Mach  number  for  group  (a)  points  exhibited  a  much  larger  variation  :  between  30  and 
80,  with  the  higher  Mach  numbers  obtained  at  further  target  points. 

These  results  are  interpreted  as  follows.  The  ambient  no-collision  probability  exp[  -  r^CS)!  IS 
sufficiently  close  to  unity,  so  that  in  an  order-of-magnitude  analysis  such  as  the  present  one.  we  may 
disregard  this  factor.  If  the  velocity  ratio  in  the  *\ik(S)  integral  of  Eq.  (3-1)  is  assumed  to  be  unity  (its 
average  value  for  group  (a)  points  is  about  1.4),  then  the  differential  in  the  flux  S  integration  becomes 
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(TN(S)dS  =  drjik(S).  This  implies  that  the  flux  S  integration  results  in  some  average  value  of  the  only 
remaining  factor  :  hjPik(S,  —  ft).  Since  the  £1  integration  introduces  a  factor  of  order  unity,  the 
order-of-magnitude  estimate  for  the  arriving-to-incident  flux  ratio  is  [hjP|k(S,  —  n)jav  .  The  value  of 
this  estimate  is  [hjPik(S,—  £2)lav=  .226  *  .09*. 02,  which  is  about  equal  to  the  actual  flux  ratio  for 
target  points  in  group  (a). 

When  an  exhaust  flow  and  orbital  parameters  (velocity  and  attitude)  are  specified,  Pik(S,  —  H) 
depends  on  the  choice  of  molecular  collision  model  (we  chose  hard  spheres),  while  hj  is  uniformly 
constant.  The  foregoing  reasoning  thus  establishes  the  collision  model  as  a  significant  factor  in 
determining  ambient  scattering  flux  levels,  to  the  extent  that  Pjk(S,  —  ft)  is  sensitive  to  the  choice  of 
model. 


Table  4-1.    Typical  Operating  Conditions  of  HF/DF  Laser  Exhaust 

Mole  fractions        [H  ]  =  .091        [HF]  =  .091        [H2]  =  .104        [DF]  =  .135        [He]  =  .579 

Average  molecular  weight  7.14 

Specific  heats  ratio  1.54 

Stagnation  temperature  and  density  1400  (K) 

.0075  (kg/m3) 

Exit  Mach  number  4.0 

Molecular  diameter  (hard  spheres)  2.5x10       (m) 

Spacecraft  diameter  5.0  (m) 

Nozzle  aperture  1.0   (m) 
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Figure  4-1.      Variation  of  Return  Flux  with  Target  Point  (Xs).    Target  Point  at  Incidence-Plane 
(<PA=0)  and  Constant  Incidence-Angle  (\\i  .=20°). 
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Figure  4-2.      Variation  of  Return  Flux  with  Ambient  Incidence  Angle  (yA).    Fixed  Target  Point 
(X  =  1  m)  Located  at  Incidence-Plane  ((pA=0). 
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Figure  4-3.      Variation  of  Return  Flux  with  Ambient  Azimuth  Angle  (cpA).     Fixed  Target  Point 
(Xs=  I  m)  and  Ambient  Incidence  Angle  (\j/A  =  20°). 
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5.       SPACECRAFT  CHARGING 

Spacecraft  charging  is  a  major  concern  to  spacecraft  designers,  particularly  for  missions  in  GEO 
and  to  a  lesser  extent  also  in  LEO.  The  exhaust  plume  of  an  HF/DF  laser  operating  in  the 
ionosphere  (300  to  1000  km  altitude)  may  modify  significantly  the  pre-firing  charging  pattern  of  the 
spacecraft.   Two  classes  of  effects  may  lead  to  charging  modification;  they  are: 

(a)  The  exhaust  contains  large  concentrations  of  HF  and  DF  molecules  which  are  highly 
electronegative.  They  may  be  readily  ionized  by  environmental  electrons  and  change  the  existing 
spacecraft  charging  pattern. 

(b)  When  the  spacecraft  is  oriented  obliquely  relative  to  its  orbital  velocity  and  the  ambient  plasma 
impinges  at  the  plume  boundary,  the  plume  will  cast  a  "shadow"  on  the  downstream  side, 
leading  to  a  very  dissimilar  charging  fluxes  on  the  upstream  and  downstream  halves  of  the 
spacecraft. 

The  knowledge  gained  in  analyzing  the  ambient  scattering  effect  can  be  applied  to  the  assessment 
of  the  effects  of  ionospheric  plasma  on  spacecraft  charging.  We  first  consider  the  upstream  side  of 
the  spacecraft  as  mentioned  in  (a)  above. 

We  contend  that  the  exhaust-plasma  interaction  will  not  drastically  alter  the  charging  pattern  of 
the  upstream  half.  This  assessment  is  established  as  follows.  Consider  the  fact  that  ionospheric 
plasma  has  a  particle  number  density  no  higher  than  1012  (m"3)  and  energy  per  particle  of  at  most 
1  (ev)  (excluding  the  auroral  plasma  of  polar  zones  or  events  of  sun  storms,  where  the  energy  per 
particle  is  much  hisher).  Significantlv,  the  Debve  length  at  the  highest  plasma  densitv  is  verv  small  : 
only  about  10°  (m);  the  largest  Debye  length  in  the  ionosphere  is  10  (m)  [9]  .  Ion  thermal  velocity 
is  typically  lower  than  orbital  velocity,  but  electron  velocity  is  considerably  higher  than  orbital 
velocity  (at  1  ev  the  electron  velocity  is  about  U  =6x  10s  m/sec).  Hence,  ions  would  typically 
impinge  at  the  plume  as  a  uniform  ion  beam  with  the  orbital  velocity  (like  ambient  molecules),  while 
electrons  are  expected  to  impinge  at  the  plume  with  their  random-oriented  thermal  velocity. 

In  view  of  the  results  of  ambient  scattering  analysis  (Ch.  3  and  4  above),  and  since  ions  are  subject 
to  similar  collision  process  with  exhaust  molecules  as  neutrals,  ions  will  be  stopped  at  the  plume 
fringes  much  like  ambient  molecules.  By  virtue  of  the  small  Debye  length  (typically  much  smaller 
than  the  stopping  distance),  electrons  would  not  penetrate  any  further  than  ions,  regardless  of  their 
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collision  cross-section  with  exhaust  molecules.  The  familiar  plasma  sheath  that  forms  at  a  solid 
surface,  is  hence  replaced  at  the  plume/plasma  boundary  by  a  typically  neutral  layer  whose  thickness 
is  of  the  order  of  an  ion/neutral  mean  free  path,  but  much  larger  than  the  Debye  length.  Only  at  the 
upper  altitude  range  of  the  ionosphere  does  the  Debye  length  become  comparable  to  a  plume 
boundary  mean  free  path  (about  .1  m),  but  there  plasma  density  and  flux  are  several  orders  of 
magnitude  lower  and  charging  modification  is  not  likely  to  be  significant  at  the  relatively  short  firing 
duration  of  about  5  minutes. 

Elastically  scattered  ions  can  be  deflected  towards  the  spacecraft  as  a  result  of  elastic  collisions 
with  exhaust  molecules,  much  like  neutrals.  Referring  to  our  analysis  of  the  return-to-ambient  flux 
ratio  (Ch.  4  above),  it  is  clear  that  the  relevant  ratio  here  will  be  about  1/4tt,  i.e.,  of  the  order  of  10% 
(this  is  due  primarily  to  the  role  played  by  the  elastic  directional  distribution  function  —  see  Ch.  4). 
A  change  in  the  plasma-to-surface  current  of  that  order  is  hence  possible,  but  unlikely  to  affect 
spacecraft  design  or  operation  significantly.  The  reason  is  that  a  design  capable  of  smoothing  away 
the  inhomogeneous  charge  flux  at  oblique  attitudes,  will  not  be  sensitive  to  a  change  in  flux  pattern 
of  the  order  of  10%  (in  other  words,  potential  differences  may  be  amplified  by  10%,  which  is  hardly 
likely  in  a  sound  design  to  bring  about  arcing  or  other  threshold  phenomena). 

Another  effect  which  may  potentially  be  significant  in  the  upstream  half  is  generation  of 
electronegative  species  (HF~ ,  DF~ )  by  plasma  electrons  impinging  at  the  plume.  In  the  sequel,  we 
examine  the  magnitude  of  this  effect,  concluding  that  it  is  negligible. 

This  estimate  is  best  done  by  considering  N     ,  which  is  the  rate  of  production  of  HF      and  DF 
per  unit  volume,  at  a  typical  point  in  the  exhaust  where  local  Mach  number  is  M  =  30  (this  is 
typically  the  lowest  average  Mach  number  for  the  plume  region  where  ambient  scattering  takes  place 
-    see  Ch.  4  above).    Since  energy  is  released  by  the  electronegative  ion  formation,  the  reaction 
involves  a  third  body  as  follows  : 

HF  +  e~   +  M     ->     HF"   +  M  (5-1) 

where  M  is  the  third  body  molecule.  We  assume  a  simplified  classical  kinetic  model  for  this  reaction, 
as  follows.  The  pair  HF/M  collide  with  a  frequency  proportional  to  the  local  number  density  and  HF 
molar  fraction,  and  to  the  average  relative  velocity.  An  electronegative  ion  formation  can  occur  only 
if  an  electron  collides  with  the  pair  during  their  collision,  which  lasts  t  =  D/Cr,  where  Cr  is  the 
average  relative  pair  velocity.    Based  on  this  classical  model,  and  assuming  the  same  cross-section  for 
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electronegative  ion  formation  as  for  elastic  HF/M  collisions,  the  volume  rate  of  electronegative  ion 
generation  is  given  by  : 

N  "  =  (7iD3N)  Nh  (7tD2UeNe)  (5-2) 

where  (7iD3N)  is  the  probability  that  a  certain  HF  or  DF  molecule  will  be  in  contact  (D  being 

molecular  hard-sphere  diameter)  with  any  other  exhaust  molecule  (whose  number  density  is  N). 

When  (7iD3N)  is  multiplied  by  hN,  where  h  is  the  HF  +  DF  molar  fraction  (Table  4-1),  the  combined 

term  reads  as  the  number  of  colliding  HF/M  pairs  per  unit  volume.    Assuming  the  electronegative 

formation  cross-section  is  also  7tD2,  the  factor  ttD2U  N    where  U    and  N    are  electron  velocitv  and 

e    e  e  e 

number  density,  renders  the  expression  for  electronegative  generation  rate  per  unit  volume.  We  note 
that  Cf  cancels  out  in  deriving  Eq.  (5-2),  so  that  N  does  not  depend  on  temperature.  This  supports 
the  use  of  the  kinetic  approximation  in  regions  of  continuum  breakdown  (plume  fringes  are  such 
regions). 

• » « 

How  is  the  relative  magnitude  of  N      decided?   To  do  that  we  multiply  N      by  X=  l/rcDzN,  which 

is  the  mean  free  path  for  a  fast  moving  particle  that  penetrates  the  plume.   This  expression  is  justified 

by  the  fact  that  most  incident  particles  do  collide  within  a  distance  of  order  X,  and  when  the  particles 

are  plasma  ions,  electrons  will  adhere  to  ion  spatial  distribution  by  virtue  of  the  small  Debye  length 

(smaller  than  X).    Thus,  XN       is  the  rate  of  electronegative  ion  generation  per  unit  area  of  plume 

boundary.   The  ratio  P       between  this  rate  and  the  incident  electron  flux  is  : 

p"   =  >.N"/NeUe  =  (7tD3N)h  =  2.2  x  10"10  (5-3) 

where  N  =  2  x  1019  (m"3)  which  corresponds  to  Mach  number  M  =  30  in  the  typical  case  (Table  4-1). 
The  fraction  of  electron  flux  captured  by  HF  and  DF  exhaust  molecules  to  form  electronegative  ions 
is  so  small  (due  to  the  pair-formation  term  (7tD3N)),  that  it  cannot  appreciably  alter  the  charging  flux 
distribution  at  the  spacecraft  surface. 

Another  possible  effect  is  the  recoil  of  HF  "    or  DF      that  occurs  due  to  energy  released  in  the 

electronegative  formation  reaction.    The  recoiling  species  might  conceivably  reach  the  surface  and 

contaminate   it.     The   magnitude   of  the   recoil   flux  is   certainly  no   larger  than   p~UN  =1300 

-2-1  _ 

(m    sec   ),  where  we  assume  the  worst  case  flux  :    Ne=  101",  Ue  =  6  *  10s  (m/sec)  which  corresponds 

to  about  1  ev  energy  per  electron.    This  flux  level  is  about  3  x  10      monolayers  of  HF      and  DF 

per  hour,  so  that  its  contribution  to  surface  contamination  is  utterly  negligible. 
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The  second  kind  of  charging  effects  (item  (b)  above)  is  due  to  the  fact  that  the  exhaust  plume  is 
impenetrable  to  ambient  plasma  (within  a  range  of  sufficiently  small  distance  from  the  spacecraft,  so 
that  no  extensive  diluting  of  the  plume  has  taken  place).  The  downstream  half  of  the  spacecraft  in 
oblique  attitude  will  be  in  the  "shadow"  with  respect  to  incident  plasma.  As  a  first  approximation  we 
may  assume  zero  plasma  flux  at  the  shadowed  surface.  More  accurately,  this  portion  of  the 
spacecraft  will  be  subject  to  a  plasma  wake  flow.  However,  it  is  quite  difficult  to  determine  the 
charging  phenomena  that  take  place  in  such  a  wake,  as  indicated  by  a  recent  work  on  solar  sails  in 
LEO  [9]  .   Thus,  a  zero  flux  at  the  downstream  half  seems  a  practical  design  assumption. 

Can  adverse  charging  effects  occur  as  a  result  of  shadowing  the  downstream  half?  This  question 
can  be  discussed  only  qualitatively.  The  reason  is  that  a  quantitative  analysis  requires  a  lumped- 
circuit  model  of  the  spacecraft  external  surface  [10]  .  Since  such  a  concrete  design  is  not  available,  we 
can  only  discuss  this  question  qualitatively.  Obviously,  assuming  zero  flux  to  the  downstream  half 
during  the  envisioned  5  minutes  of  laser  firing  duration,  and  requiring  that  no  appreciable  voltages 
between  the  two  halves  will  evolve,  leads  to  the  stipulation  that  the  equivalent-circuit 
Capacitance  x  Resistance  should  be  much  smaller  than  the  firing  duration. 
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6.        CONCLUDING   REMARKS 

Our  major  quantitative  conclusion  is  that  for  the  relatively  high  ambient  density  assumed 
(NA=lxiOi6  molecules/m3  which  represents  Sunspot  Maximum  at  about  200  km)  and  for  the 
typical  HF/DF  laser  exhaust  (Table  4-1),  the  HF  +  DF  flux  backscattered  by  ambient  molecules  is 
several  hundred  monolayers  per  hour.  This  flux  level  may  seem  as  not  outright  negligible.  However, 
since  ambient  scattering  flux  is  proportional  to  ambient  density,  it  will  be  scaled  down  considerably  at 
the  lower  ambient  densities  of  higher  orbital  altitudes. 

The  operational  scenario  for  HF/DF  laser  envisions  4  or  5  minutes  total  operating  time;  hence  the 
contamination  by  ambient  scattering  may  not  be  serious  due  to  short  operating  time. 

The  effects  of  laser  exhaust  plume  on  spacecraft  charging  in  the  ionosphere  were  examined.  It  was 
concluded  that  the  rate  of  electronegative  (HF  ~  and  DF  ~ )  production  by  impinging  electrons  was 
negligible;  the  low  rate  is  a  consequence  of  the  assumption  that  a  third  body  is  required  to  interact 
simultaneously  with  the  HF/e  or  DF/e  pair.  No  significant  modification  of  charging  pattern  is 
anticipated.  However,  at  oblique  orbital  attitudes,  the  downstream  half  of  the  spacecraft  will  be 
shadowed  from  the  oncoming  ambient  plasma.  This  fact  has  to  be  reckoned  with  in  designing  a  ring- 
symmetric  laser  spacecraft. 

The  emphasis  in  this  work  was  on  the  method  rather  than  on  results.  The  first-collision  model  was 
demonstrated  to  be  simple  to  implement  in  a  code.  It  is  considerably  simpler  than  the  more  general 
and  potentially  more  accurate  Monte  Carlo  methods  commonly  used  for  simulating  rarefied  flows  [8]  . 
We  found  out  that  the  molecular  collision  model  was  all  important  in  determining  the  return  flux 
level,  which  is  hardly  surprising  for  scattering  by  single  collision.  For  the  same  reason,  the  collision 
model  would  also  be  dominant  in  a  Monte  Carlo  simulation  of  the  ambient  scattering  process. 

If  and  when  a  mathematical  accuracy  of  the  first-collision  approximation  is  established  for  hard- 
spheres,  it  might  be  possible  to  determine  a  realistic  collision  model  by  comparing  computed  results 
with  measurements. 

This  accuracy  may  be  established  in  either  of  two  ways.  One  way  is  by  comparison  with  accurate 
Monte  Carlo  computations  (using  hard-spheres  collision  model).  The  other  is  to  seek  an  estimate  o[ 
the  error  incurred  by  considering  just  first  collisions  and  ignoring  all  subsequent  ones.  This  might  be 
achieved  by  accounting  for  second  collisions  in  an  extended  first-collision  model,  provided  a  simplified 
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scheme  that  will  obviate  the  need  for  increase  in  the  dimensions  of  the  numerical  flux  integration  can 
be  devised.   We  are  currently  considering  such  second-collision  approaches. 
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APPENDIX  A.      DESCRIPTION   OF  AMB  CODE 


A.l      Description  of  Subroutines 

We  provide  a  list  of  the  subroutines  in  the  ambient  flux  integration  code  AMB  for  ring- symmetric 
cylidrical  spacecrafts.  Each  subroutine  is  briefly  described.  Statements  are  identified  by  the 
FORTRAN  statement  number  (columns  1  through  5). 

MAIN  PROGRAM  The  300  loop  is  intended  to  enable  several  (NCASE)  reruns  with  various  data 
in  each,  all  in  a  single  run.  Upon  calling  INIDAT1,  parameters  depending  on  data  defined 
in  INIDAT  are  re-computed.  The  200  loop  is  over  various  XSV(NX)  target  points.  In  the 
20  loop  the  flux  integration  begins  :  FLUXC  is  for  particle  flux  and  FLXU2C  is  for  the 
rms  of  velocity  of  return  flux  molecules.  All  the  MAX  suffixed  parameters  denote  values  at 
which  the  integrand  had  the  largest  value. 

The  actual  flux  integration  commences  at  statement  1  for  the  summation  over  strips  of 
constant  RF.  This  summation  is  terminated  when  convergence  is  attained  (to  within 
EPSR).  The  inner  loop  2  is  over  azimuth  angle  PHI.  Note  that  the  target  points  are 
generally  not  in  the  plane  of  incidence  (PHIA.NE.0),  so  that  no  symmetry  can  be  assumed 
in  the  PHI  integration,  and  it  is  performed  twice  in  order  to  cover  the  entire  range  in  PHI 
(IPAR=  1  for  PHI.GT.0,  IPAR=2  for.PHI.LT.0).  The  flux  integration  along  the  line-of- 
sight  is  done  by  calling  FLUX. 

INIDAT  Initialization  of  data.  There  is  no  input  file  for  this  code.  INIDAT1  is  for  parameters 
computed  from  the  data  defined  by  calling  INIDAT. 

SOF  Stopping  routine,  called  when  an  error  is  detected.    Here  we  also  trigger  a  system  error  by 

computing  DSQRT(-l),  in  order  to  obtain  a  calling  sequence  printout  by  the  operating 
system. 

FLUX        This  routine  calls  SUMT  for  flux  integration  of  one  exhaust  species  at  a  time. 

LIMIT  Here  we  compute  the  point  of  intersection  of  the  line-of-sight  with  the  leading 
characteristic  cone.  If  they  do  not  intersect,  the  distance  of  the  intersecting  point  TLIM  is 
set  to  a  verv  larse  number. 
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SUMT  This  is  the  line-of-sight  integration  routine.  Runge-Kutta  scheme  is  used  (even  though  an 
explicit  integral  is  computed).  Note  that  ETAK  and  ETAIK  have  to  be  computed  through 
a  separate  integration  at  each  point  of  the  line-of-sight  integration.  The  integration  step 
DT  (T=S/RF)  is  re-adjusted  at  each  integration  step.  The  integration  is  terminated  when 
convergence  is  attained  (to  within  EPST). 

FETA  Here  the  integrand  for  the  line-of-sight  flux  integration  is  evaluated.  The  hard-spheres 
collision  model  is  used  to  determine  the  post-collision  directional  distribution  factor  PIK. 
The  flux-average  of  any  variable  (such  as  UIK**2  in  present  version),  can  be  computed  by 
summing  it  multiplied  by  flux  and  subsequently  dividing  by  the  total  arriving  flux  (see  loop 
31  in  MAIN  PROGRAM). 

PATH1K  Here  the  molecular  thickness  ETAIK  of  the  I  exhaust  species  scattered  by  the  K  ambient 
molecule,  is  computed  by  integration  along  the  line-of-sight. 

FT      This  routine  computes  the  integrand  for  the  ETAIK  integration  in  PATHIK. 

PATHK  The  analog  to  PATHIK  for  K  ambient  molecule.  TAU  is  the  normalized  integration 
variable  along  the  trajectory  of  the  penetrating  ambient  molecule.  Note  that 
SHADOW  =  .TRUE,  when  the  trajectory  passes  through  the  cylindrical  spacecraft  surface 
before  entering  the  fan. 

FTAU        Computes  the  integrand  for  the  ETAK  integration  in  PATHK. 

FAN  Computes  the  fan  coordinates  PSI,  XP,  YP  for  a  point  on  the  line-of-sight.  It  is  used  to 
determine  the  Mach  number  and  flow  angle  from  the  power-law  approximation  (see 
MATCH). 

FANT        Computes  the  fan  coordinates  PSI,  XP,  YP  for  a  point  on  the  ambient  molecule  trajectory. 

HMSET  Prepares  the  vector  HMV(I)  which  is  the  value  of  the  H(M)  integral  at  a  set  of  Vlach 
number  values  (equally  spaced  in  inverse  Mach  number).  This  vector  is  used  to  compute 
H(M)  for  an  arbitrary  M  (see  H INTER),  since  this  function  is  needed  in  the  power-law 
approximation  of  flow  in  a  ring-symmetric  fan.  Subsequent  routines  MFUNC,  HINTER, 
MATCH  and  AREAF  are  all  used  to  implement  this  approximation. 
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MFUNC      Computes  the  integrand  for  the  H(M)  integration  in  HMSET. 

HINTER  Computes  H(M)  for  a  given  M,  from  HMV(I)  by  linear  interpolation.  Note  that  the 
interpolation  is  done  with  inverse  Mach  number  as  the  independent  variable. 

MATCH  Here  the  approximation  to  the  "inverse  problem"  of  finding  the  Mach  number  at  a  single 
point  in  the  ring-fan  is  implemented.  An  iteration  scheme  is  used  to  determine  the  fan 
characteristic  passing  through  the  given  point  [2]  . 

AREAF  iMach  number  is  computed  from  value  of  area  ratio  function.  Newton-Raphson  iterations 
are  used. 
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A.2      Listing  of  AMB  code 

C$0PTI0NS  LIST  AMB0001 

C   AMBIENT.   SCATTERING  FROM  A  RING  PLUME  BY  AMBIENT  AIR.  AMB0002 

IMPLICIT  REAL*8(A-H,0-Z)  AMB0003 

COMMON  /GAMA/G,G1,G2,G3, G4, G5, G6 ,G7 ,G8,G9, G10, Gil , G12,G13, G14, G15, AMB 00 04 

1             G16,G17,G18,G19,G20  AMB0005 

COMMON  /PAR/CO, ENO, EMI, D, SIGMA, TLIM,DR0, EL 0,Q0, TO, FACT, ALOGF,  AMB0006 

1            DPSIO,DTMAX,DETAO,ETALIM,XSI,XSF  AMB0007 

COMMON  /NPAR/NPHI,IPAR,NP,NR,NX,NXS,NS,NSPEC,NS1,NS2,NTAU0,NETA0,  AMB0008 

1             NAMB,NCASE,ICASE,IFAN  AMB0009 

COMMON  /GE0M/APF,PAI,PAI2,W,SW,CW,BETA,SBETA,CBETA,PSI1,SPSI1,     AMB0010 

1  CPSI1,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,A0,RF,XF,YF,ZF,AMB0011 

2  PHISOF,PHIF,SPHIF,CPHIF,DYMIN,RMIN,XS,DIST,X0,Y0,Z0,  AMB0012 

3  DYO,DEG,PSIN,ST1,CT1,OMEGX,OMEGY,OMEGZ,XSV(21)  AMB0013 
COMMON  /EPSIL/EPSETA,EPST,EPSR  AMB0014 
COMMON  /EXTREM/TEXT(5),ETAEXT(5),ETAKXT(5),PHIEXT(5),  AMB0015 

1  PSIEXT(5),EMEXT(5),FEXT(5),WEXT(5),  AMB0016 

2  TMAX(5),ETAKMX(5),ETAMAX(5),PSIMAX(5),  AMB0017 

3  EMMAX(5),FMAX(5),  AMB0018 

4  RFMAXC5),PHIFMX(5),PHIMAX(5),WMAX(5)  AMB0019 
COMMON  /COUNTS/ICONTC, ICONTT, ICNTOT, ICNTMX, IQT0TC5) , ISHADC5)  AMB0020 
COMMON  /SPEC/WAV,XCC5),WC(5),WRC(5),XNAME(5),QFC(5),QDC(5),  AMB0021 

1             QU2C(5),FLUXC(5),0MEGA(5),FLXU2C(5),URMSC(5)  AMB0022 

COMMON  /AMBIEN/ENA,UA,PSIA,PHIA,HAC3),WA(3),  AMB0023 

1               UAX,UAY,UAZ,AA,BA,CA,RA,XA,YA,ZA, SHADOW  AMB0024 

COMMON  /POINT/XP,YP,XCOR,YCOR  AMB0025 

LOGICAL  SHADOW  AMB0026 

DIMENSION  DSUMF(5),DSUMD(5),DSUMAX(5),DSUMU2(5)  AMB0027 

NCASE=1  AMB0028 

DO  300  ICASE=1,NCASE  AMB0029 

CALL  INIDAT  AMB0030 

GO  TO  (301,302,303,304,305,306,307,308,309,310,  AMB0031 

1       311, 312, 313, 314, 315, 316, 317, 318, 319, 320), ICASE  AMB0032 

301  CONTINUE  AMB0033 
IFAN=1  AMB0034 
NXS=3  AMB0035 
XSI=0.1D0  AMB0036 
GO  TO  399  AMB0037 

302  CONTINUE  AMB0038 
PHIA=20.D0/DEG  AMB0039 
GO  TO  399  AMB0040 

303  CONTINUE  AMB0041 
PHIA=50.D0/DEG  AMB0042 
GO  TO  399  AMB0043 

304  CONTINUE  AMB0044 
PHIA=75.D0/DEG  AMB0045 
GO  TO  399  AMB0046 

305  CONTINUE  AMB0047 
PHIA=100.D0/DEG  AMB0048 
GO  TO  399  AMB0049 

306  CONTINUE  AMB0050 
PHIA=125.D0/DEG  AMB0051 
GO  TO  399  AMB0052 

307  CONTINUE  AMB0053 
PHIA=150.D0/DEG  AMB0054 
GO  TO  399  AMB0055 

308  CONTINUE  AMB0056 
PHIA=175.D0/DEG  AMB0057 
GO  TO  399  AMB0058 

309  CONTINUE  AMB0059 
GO  TO  399  AMB0060 

310  CONTINUE  AMB0061 
GO  TO  399  AMB0062 

311  CONTINUE  AMB0063 
GO  TO  399  AMB0064 

312  CONTINUE  AMB0065 
GO  TO  399  AMB0066 

313  CONTINUE  AMB0067 
GO  TO  399  AMB0068 

314  CONTINUE  AMB0069 
GO  TO  399  AMB0070 

315  CONTINUE  AMB0071 
GO  TO  399  AMB0072 

34 


SCRIPT    Al 

316 

CONTINUE 

GO  TO  399 

317 

CONTINUE 

GO  TO  399 

318 

CONTINUE 

GO  TO  399 

319 

CONTINUE 

GO  TO  399 

320 

CONTINUE 

GO  TO  399 

399 

CONTINUE 

PRINT  101 

101 

FORMATCl') 

C 

AMB0073 
AMB0074 
AMB0075 
AMB0076 
AMB0077 
AMB0078 
AMB0079 
AMB0080 
AMB0081 
AMB0082 
AMB0083 
AMB0084 
AMB0085 
AMB0086 

CALL  INDAT1  AMB0037 

C  AMB0088 

DO  200  NX=1,NXS  AMB0089 

XS=XSV(NX)  AMB0090 

C   CX0,Y0,Z0)  IS  THE  POINT  AT  WHICH  FLUX  AND  DENSITY  ARE  COMPUTED.       AMB0091 
C   THE  NORMAL  TO  THE  SURFACE  AT  (X0,Y0,Z0)  IS  PARALLEL  TO  Y-AXIS.        AMB0092 
X0=XS  AMB0093 

Y0=A0  AMB0094 

Z0=0.  AMB0095 

DO  20  NS=NS1,NS2  AMB0096 

FLUXC(NS)=0.  AMB0097 

FLXU2C(NS)=0.  AMB0098 

OMEGA(NS)=0.  AMB0099 

DSUMAX(NS)=0.  AMB0100 

IQTOT(NS)=0  AMB0101 

ISHAD(NS)=0.  AMB0102 

TMAX(NS)=-1.D  44  AMB0103 

ETAKMX(NS)=-1.D  44  AMB0104 

PHIMAX(NS)=-1.D  44  AMB0105 

PHIFMX(NS)=-1.D  44  AMB0106 

WMAX(NS)=-1.D  44  AMB0107 

PSIMAX(NS)=-1.D  44  AMB0108 

ETAMAX(NS)=-1.D  44  AMB0109 

RFMAX(NS)=-1.D  44  AMB0110 

EMMAX(NS)=-1.D  44  AMB0111 

FMAX(NS)=-1.D  44  AMB0112 

20  CONTINUE  AMB0113 
RN=RMIN  AMB0114 

APF=A0-0.5D0*DR0*SPSIF  AMB0115 

NR=0  AMB0116 

1     NR=NR+1  AMB0117 

DR=DR0  AMB0118 

DR=DR0*(APF/A0)  AMB0119 

C     DR=DR0*(1.D0+0.4D0*DR0/XS)**NR  AMB0120 

RF=RN+DR/2.D0  AMB0121 

APF=AO+RF*SPSIF  AMB0122 

PHISOF=DACOS(A0/APF)  AMB0123 

C     DPHIO  =  0.1DO  AMB0124 

C     NPHI=PHISOF/DPHI0+2  AMB0125 

DPHI=PHISOF/DBLE(NPHI)  AMB0126 

DO  21  NS=NS1,NS2  AMB0127 

DSUMF(NS)=0.  AMB0128 

DSUMU2(NS)=0.  AMB0129 

DSUMD(NS)=0.  AMB0130 

21  CONTINUE  AMB0131 
DOMEGR=0.  AMB0132 
DO  2  NP=1,NPHI  AMB0133 
DO  2  IPAR=1,2  AMB0134 
PHIF=(DBLE(NP)-0.5D0)*DPHI  AMB0135 
IFCIPAR.EQ.2)  PHIF=-PHIF  AMB0136 

C  AMB0137 

CALL  FLUX  AMB0138 

C  AMB0139 

CROSS1=OMEGY  AMB0140 

CR0SS2=(SPSIF)*(-0MEGX)+(-CPSIF*CPHIF)*(-0MEGY)+  AMB0141 

1        (-CPSIF*SPHIF)*(-OMEGZ)  AMB0142 

IF(CROSS1.LE.O.)  AMB0143 

1CALL  SOFCDIRECTION  COSINE  OF  SURFACE  NORMAL  SHOULD  BE  POSITIVE')  AMB0144 
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IFCCROSS2.LE.0.) 
1CALL    SOFCNORMAL    TO    LIMITING   CONE   HAS   NEGATIVE   PROJECTION    ON 
10F-SIGHT') 

D0MEGA=CR0SS2*DPHIXAPF*DR/DISTXX2 

DOMEGR=DOMEGR+DOMEGA 

DO    24    NS=NS1,NS2 

DSUMF(NS)=DSUMFCNS)+DOMEGAXQFC(NS)*CROSSl 

DSUMU2(NS)=DSUMU2(NS)+D0MEGA*QU2C(NS)*CR0SS1 

IF(DSUMAXCNS) .GT.DOMEGA*QFC(NS)*CROSSl)  GO  TO  24 

DSUMAX(NS)=DOMEGA*QFC(NS)*CROSSl 

TMAX(NS)=TEXT(NS) 

ETAKMX(NS)=ETAKXTCNS) 

PHIMAX(NS)=PHIEXTCNS)*DEG 

PHIFMX(NS)=PHIF*DEG 

WMAX(NS)=WEXT(NS)*DEG 

PSIMAX(NS)=PSIEXT(NS)*DEG 

ETAMAX(NS)=ETAEXT(NS) 

RFMAX(NS)=RF 

EMMAX(NS)=EMEXT(NS) 

FMAX(NS)=QFC(NS)*XC(NS)XQO 

CONTINUE 

CONTINUE 

DO  26  NS=NS1,NS2 

FLUXC(NS)=FLUXC(NS)+DSUMF(NS) 

FLXU2C(NS)=FLXU2C(NS)+DSUMU2(NS) 

OMEGA(NS)=OMEGA(NS)+DOMEGR 

CONTINUE 

RN=RN+DR 

IFCNR.LE.2)  GO  TO  1 

IFCNR.GT.99)  GO  TO  10 

DO  27  NS=NS1,NS2 

IF(FLUXC(NS).EQ.O.)  GO  TO  27 

ERR=(DSUMF(NS)/FLUXC(NS))/DOMEGR 

IF(ERR.GT.EPSR)  GO  TO  28 

CONTINUE 

GO  TO  10 

CONTINUE 

GO  TO  1 

CONTINUE 

DO  31  NS=NS1,NS2 

FLUXC(NS)=XC(NS)*FLUXC(NS)*QO 

OMEGA(NS)=OMEGA(NS)/(2.D0*PAI*DCOS(PSIF/2.D0)**2) 

FLXU2C(NS)=XC(NS)*FLXU2C(NS)*Q0 

URMSC(NS)=0. 

IFCFLUXC(NS) .EQ.O.)  GO  TO  31 

URMSC(NS)=DSQRT(FLXU2C(NS)/FLUXC(NS)) 
:   AVERAGE  EM   (SEE  FETA) 

URMSC(NS)=       FLXU2C(NS)/FLUXC(NS) 
31    CONTINUE 

PRINT  ll,NX,NR,XS,RF,DR,PHISOF*DEG 

11  FORMATC///1X, 'NX,NR,XS,RF,DR,PHIS0F=',2I4,3D13.4,F8.4, 
1        3X,'FLUX  AND  EXTREMA  VALUES,  ALL  SPECIES:'/) 

PRINT  12 

12  FORMATC/1X,1  NAME  ','  IQTOT',1  ISHAD', 


24 
2 


26 


27 

28 


10 


NAME  ', ' 

FMAX    ' 

ETAKMX' 

EMMAX' 

URMSC 

NS2 


13 
14 
200 


FORMATC/1X, • 

1  • 

2  • 

3  ' 

4  » 
DO  14  NS=NS1 
DLF=0. 
IFCFLUXC(NS) .NE.O) 

1DLF=DLOG10(FLUXC(NS))+100 

IDLF=DLF 

DLF=DLF-DBLE(IDLF) 

PRINT  13,XNAME(NS),IQT0T(NS) 
1 
2 
3 
4 


OMEGA* 
ETAMAX' 

RFMAX' 
FLUXC 


TMAX' 
PSIMAX' 
PI-WMAX' 
LOG'/) 


F0RMAT(1X,A6,2I6 

CONTINUE 

CONTINUE 


D0+1.D-11 


ISHAD( NS),FMAX(NS), OMEGA  (NS), 
TMAXC  NS) , ETAKMX( NS ) , ETAMAXC  NS) 
PS IMAX(NS), EMMAX ( NS) , RFMAX( NS) 
180.D0-WMAX(NS),URMSC(NS), 
FLUXC(NS),DLF 


D10.3,4F8.4,4F8.1,F8.2,D10.3,  V,F4.2) 


AMB0145 
LINE-AMB0146 
AMB0147 
AMB0148 
AMB0149 
AMB0150 
AMB0151 
AMB0152 
AMB0153 
AMB0154 
AMB0155 
AMB0156 
AMB0157 
AMB0158 
AMB0159 
AMB0160 
AMB0161 
AMB0162 
AMB0163 
AMB0164 
AMB0165 
AMB0166 
AMB0167 
AMB0168 
AMB0169 
AMB0170 
AMB0171 
AMB0172 
AMB0173 
AMB0174 
AMB0175 
AMB0176 
AMB0177 
AMB0178 
AMB0179 
AMB0180 
AMB0181 
AMB0182 
AMB0183 
AMB0184 
AMB0185 
AMB0186 
AMB0187 
AMB0188 
AMB0189 
AMB0190 
AMB0191 
AMB0192 
AMB0193 
AMB0194 
AMB0195 
AMB0196 
AMB0197 
AMB0198 
AMB0199 
AMB0200 
AMB0201 
AMB0202 
AMB0203 
AMB0204 
AMB0205 
AMB0206 
AMB0207 
AMB0208 
AMB0209 
AMB0210 
AMB0211 
AMB0212 
AMB0213 
AMB0214 
AMB0215 
AMB0216 
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PRINT  102  AMB0217 

102   F0RMAT(///1X,'END   RING   RUN',///)  AMB0218 

300   CONTINUE  AMB0219 

STOP  AMB0220 

END  AMB0221 

SUBROUTINE  INIDAT  AMB0222 

IMPLICIT  REAL*8(A-H,0-Z)  AMB0223 

REAL*8  LAMDA0,LAMDA1  AMB0224 

CHARACTER*8  XNAME,XNAMED  AMB0225 
COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15,AMB0226 

1             G16,G17,G18,G19,G20  AMB0227 

COMMON  /PAR/CO, ENO, EMI, D, SIGMA, TLIM,DRO, EL 0,QO, TO, FACT, ALOGF,  AMB0228 

1  DPSIO,DTMAX,DETAO,ETALIM,XSI,XSF  AMB0229 
COMMON  /NPAR/NPHI,IPAR,NP,NR,NX,NXS,NS,NSPEC,NS1,NS2,NTAU0,NETA0,  AMB0230 

1  NAMB,NCASE,ICASE,IFAN  AMB0231 
COMMON  /GE0M/APF,PAI,PAI2,W,SW,CW,BETA,SBETA,CBETA,PSI1,SPSI1,     AMB0232 

1  CPSI1,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,A0,RF,XF,YF,ZF,AMB0233 

2  PHISOF,PHIF,SPHIF,CPHIF,DYMIN,RMIN,XS,DIST,X0,Y0,Z0,  AMB0234 

3  DYO,DEG,PSIN,ST1,CT1,OMEGX,OMEGY,OMEGZ,XSV(21)  AMB0235 
COMMON  /EPSIL/EPSETA,EPST,EPSR  AMB0236 
COMMON  /EXTREM/TEXT(5),ETAEXT(5),ETAKXT(5),PHIEXT(5),  AMB0237 

1  PSIEXT(5),EMEXT(5),FEXT(5),WEXT(5),  AMB0238 

2  TMAX(5),ETAKMX(5),ETAMAX(5),PSIMAX(5),  AMB0239 

3  EMMAX(5),FMAX(5),  AMB0240 

4  RFMAX(5),PHIFMX(5),PHIMAX(5),WMAX(5)  AMB0241 
COMMON  /C0UNTS/IC0NTC,IC0NTT,ICNT0T,ICNTMX,IQT0T(5),ISHAD(5)  AMB0242 
COMMON  /SPEC/WAV,XC(5),WC(5),WRCC5),XNAME(5),QFC(5),QDC(5),  AMB0243 

1             QU2C(5),FLUXC(5),0MEGA(5),FLXU2C(5),URMSC(5)  AMB0244 

COMMON  /AMBIEN/ENA,UA,PSIA,PHIA,HA(3),WA(3),  AMB0245 

1               UAX,UAY,UAZ,AA,BA,CA,RA,XA,YA,ZA, SHADOW  AMB0246 

COMMON  /POINT/XP,YP,XCOR,YCOR  AMB0247 

LOGICAL  SHADOW  AMB0248 

DIMENSION  XCD(5),WCD(5),XNAMED(5)  AMB0249 

DATA  XCD/.091D0, .091D0, .104D0, .135D0, .579D0/  AMB0250 

DATA  WCD/1.00DO,20.0DO,2.00DO,21.0DO,4.00DO/  AMB0251 

DATA  XNAMED/'    H   ',»   HF   ','   H2   • ,'   DF   ','   HE   V  AMB0252 

DATA  IFIRST/0/  AMB0253 

IFAN=2  AMB0254 

PAI=4.D0*DATAN(1.D0)  AMB0255 

PAI2=PAI/2.D0  AMB0256 

DEG=180.D0/PAI  AMB0257 

AR=8.3143D3  AMB0258 

AV=6 . 022D  26  AMB0259 

C   OMEGAC=0.5  IS  FOR  HARD  SPHERE  COLLISIONS,  AMB0260 

C   AN  AVERAGE  RECOMMENDED  VALUE  IS  ABOUT  0MEGAC=0.75  AMB0261 

OMEGAC=0.5D0  AMB0262 

NSPEC=5  AMB0263 

NS1=2  AMB0264 

NS2  =  2  AMB0265 

DO  51  NS=1,NSPEC  AMB0266 

XC(NS)=XCD(NS)  AMB0267 

WC(NS)=WCD(NS)  AMB0268 

XNAME(NS)=XNAMED(NS)  AMB0269 

51    CONTINUE  AMB0270 

C   COMBINE  HF  AND  DF  MOLE  FRACTIONS  INTO  HF  FRACTION  AMB0271 

XC(2)=XC(2)+XC(4)  AMB0272 

XC(4)=0.  AMB0273 

C  AMB0274 

A0=2.5D0  AMB0275 

EM1=4.D0  AMB0276 

RHO0=0.0075D0  AMB0277 

T0=1400.D0  AMB0278 

G=1.54D0  AMB0279 

D=2.5D-10  AMB0280 

NXS=1  AMB0281 

XSI=1.0D0  AMB0282 

XSF=10.D0  AMB0283 

C   AMBIENT  AIR  AMB0284 

ENA=1.00D  16  AMB0285 

UA=8.D  3  AMB0286 

NAMB=3  AMB0287 

WA(1)=28.D0  AMB0288 
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WA(2)=32.D0 
WA(3)=16.D0 
HA(1)=1.D0 
HA(2)=0. 
HA(3)=0. 
PSIA=20.D0/DEG 
PHIA=O.OODO/DEG 
C   INTEGRATION  PARAMETERS 
NPHI=6 
NTAU0=4 
NETA0=4 
ICNTMX=100 
RMIN=0. 
DRO=0.10DO 
DPSI0=0.20D0 
DTMAX=1.0DO 
DETA0=0.50D0 
ETALIM=10.DO 
EPST=0.5D0 
EPSR=0.3DO 
FACT=1.D  20 
RETURN 

C  COMPUTATION  OF  DATA-DEPENDENT  PARAMETERS 
Cmxxxxxxxmxxxxxxxxxxxxxx^xxxmxx^xxxxx^xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

ENTRY  INDAT1 
Cxxxx^xx^xxxxx^x^xxxxxxxxxxxxsckxx^xxxxxxxxxxxx^^x^xxxxxxxxxxxxx^xxxxxx 

ALOGF=DLOG(FACT) 

WAV=0. 

DO  52  NS=1,NSPEC 

WAV=WAV+XC(NS)*WC(NS) 

52  CONTINUE 

DO  53  NS=1,NSPEC 
WRC(NS)=WC(NS)/WAV 

53  CONTINUE 

SIGMA=PAI*D**2 

EN0=RHO0*AV/WAV 

CO=DSQRT(G*AR*TO/WAV) 

XSV(1)=XSI 

IF(NXS.EQ.l)  GO  TO  12 

DXL=(DLOG(XSF)-DLOG(XSI))/(DBLE(NXS)-l.DO) 

XLI=DLOG(XSI) 

DO  11  NX=2,NXS 

XSV(NX)=DEXP(XLI+(DBLE(NX)-1.D0)XDXL) 

CONTINUE 


11 
12 


CONTINUE 

G1=(G-1.D0)/2.D0 
G2=(G+l.D0)/(2.DO*(G-l.D0)) 
G3=G/2.D0 

G4=(G+1.D0)/(G-1.D0) 
G5=DSQRT((G+1.D0)/(G-1.D0)) 
G6=1.D0/(G-1.D0) 
G7=2.D0/(G+1.D0) 

G8=(0.5D0*(G+1 . D0)**2/(G-1 . D0))**(1 .D0/(G+1 
1    ((G+1.D0)/(G-1.D0))**((G-1.D0)/(G+1.D0)) 
G9=(G+3.DO)/(2.D0*(G-l .DO)) 
G10=(7.D0-3.D0*G)/(2.D0*(G-1 .DO)) 
G11=(5.D0-3.D0*G)/(2.D0*(G-1.D0)) 
G13=(2.D0-G)/(2.D0*(G-1.D0)) 
G14=G/(2.D0*(G-1.D0)) 
G15=(G+1.D0)/(3.D0-G) 

ZETA1=G5*DATAN(DSQRT(EM1**2-1.D0)/G5) 
AMU1=DASIN(1 .DO/ EMI) 
PSI1=PAI2+AMU1 
SPSI1=DSIN(PSI1) 
CPSI1=DC0S(PSI1) 
PSIF=PAI2+AMU1+ZETA1-G5*PAI2 
SPSIF=DSIN(PSIF) 
CPSIF=DCOS(PSIF) 
TPSIF=DTAN(PSIF) 
TETA1=PSI1-AMU1 
ST1=DSIN(TETA1) 


DO))* 


AMB0289 
AMB0290 
AMB0291 
AMB0292 
AMB0293 
AMB0294 
AMB0295 
AMB0296 
AMB0297 
AMB0298 
AMB0299 
AMB0300 
AMB0301 
AMB0302 
AMB0303 
AMB0304 
AMB0305 
AMB0306 
AMB0307 
AMB0308 
AMB0309 
AMB0310 
AMB0311 
AMB0312 
AMB0313 
AMB0314 
AMB0315 
AMB0316 
AMB0317 
AMB0318 
AMB0319 
AMB0320 
AMB0321 
AMB0322 
AMB0323 
AMB0324 
AMB0325 
AMB0326 
AMB0327 
AMB0328 
AMB0329 
AMB0330 
AMB0331 
AMB0332 
AMB0333 
AMB0334 
AMB0335 
AMB0336 
AMB0337 
AMB0338 
AMB0339 
AMB0340 
AMB0341 
AMB0342 
AMB0343 
AMB0344 
AMB0345 
AMB0346 
AMB0347 
AMB0348 
AMB0349 
AMB0350 
AMB0351 
AMB0352 
AMB0353 
AMB0354 
AMB0355 
AMB0356 
AMB0357 
AMB0358 
AMB0359 
AMB0360 
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CT1=DC0SCTETA1)  AMB0361 

QO=ENA*UA  AMB0362 

LAMDA0=1.D0/(DSQRT(2.D0)XSIGMA*EN0)  AMB0363 

LAMDA1=LAMDAO*(1.DO+G1XEM1*X2)*X(G6-OMEGAC+0.5DO)  AMB0364 

AA=DCOS(PSIA)  AMB0365 

BA=DSIN(PSIA)*DCOS(PHIA)  AMB0366 

CA=DSIN(PSIA)*DSIN(PHIA)  AMB0367 

UAX=-UA*AA  AMB0368 

UAY=-UA*BA  AMB0369 

UAZ=-UA*CA  AMB0370 

XCOR=0.  AMB0371 

YCOR=A0  AMB0372 

C  AMB0373 

PRINT  201,NSPEC,XNAME  AMB0374 

201  FORMATC/1X, 'SPECIES  DATA    NSPEC=',I3/  AMB0375 
1         IX, 'SPECIES  NAMES     MK2X,  A6 ,  2X) )  AMB0376 

PRINT  202, XC  AMB0377 

202  FORMATC  IX,  'MOLE  FRACTION  XC=  ',  1KF8  .  4,  2X)  )  AMB0378 
PRINT  203, WC  AMB0379 

203  FORMATC  1X,'M0L.  WEIGHT    WC= ' , 11 (F8 . 4, 2X) )  AMB0380 
PRINT  21,AR,AV,WAV,G,RHO0,T0,EN0,C0,D  AMB0381 

21  F0RMATC/1X, 'THERMODYNAMIC  DATA'/  AMB0382 

1  IX, 'AR,AV,WAV,GAMMA=',2X,2D14.5,2F9.3/  AMB0383 

2  1X,'RHOO,TO,ENO,CO,D=',D12.4,F8.0,D13.5,2D12.4)  AMB0384 
PRINT  22,EM1,PSI1*DEG,PSIF*DEG,  AMB0385 

1  A0,LAMDA0,LAMDA1  AMB0386 

22  FORMATC/1X, 'FLOW  AND  GEOMETRY  DATA'/  AMB0387 

1  IX, 'EM1,PSI1,PSIF=',3F9.3/  AMB0388 

2  IX, »A0,LAMDA0,LAMDA1=',F9.3,2D13.4)  AMB0389 
PRINT  23,DPSI0,DTMAX,DETA0,ETALIM,DR0,RMIN,  AMB0390 

1  EPST,EPSR,  AMB0391 

2  NPHI,NTAUO,NETAO  AMB0392 


23    F0RMAT(/1X, 

1  IX, 

2  IX, 

3  IX, 

4  IX, 


24    FORMATC/1X, 

1  IX, 

2  IX, 


INTEGRATION  DATA'/  AMB0393 

DPSIO, DTMAX, DETAO, ETALIM= ' , 4F9 . 4/  AMB0394 

DR0,RMIN,=',2D13.4/  AMB0395 

EPST,EPSR=',2D12.3/  AMB0396 

NPHI,NTAUO,NETAO=f ,316)  AMB0397 


PRINT  24,ENA,UA,PSIA*DEG,PHIA*DEG  AMB0398 


ABBREVIATED  AIR  DATA'/  AMB0399 

ENA,UA=*,2D13.4/  AMB0400 

PSIA,PHIA=',2F9.1)  AMB0401 

GO  TO  (251,252),  IFAN  AMB0402 

AMB0403 

251  CONTINUE  AMB0404 
PRINT  2510,  IFAN  AMB0405 

2510  FORMATC/1X, 'RING-FAN  APPROXIMATED  AS  PLANAR.    IFAN=',I4)  AMB0406 

GO  TO  250  AMB0407 

252  CONTINUE  AMB0408 
PRINT  2520,  IFAN  AMB0409 

2520  FORMATC/1X, 'RING-FAN  APPROXIMATED  BY  MATCHED  APPROXIMATION.',  AMB0410 

1        4X, 'IFAN=',I4)  AMB0411 

250   CONTINUE  AMB0412 

PRINT  29  AMB0413 

29    F0RMATC///1X, 'END   DATA'///)  AMB0414 

IFCIFIRST.EQ.0.AND.IFAN.EQ.2)  AMB0415 

1CALL  HMSET  AMB0416 

IFCIFAN.EQ.2)  IFIRST=IFIRST+1  AMB0417 

RETURN  AMB0418 

END  AMB0419 

C$OPTIONS  LIST  AMB0420 

SUBROUTINE  SOF(ISTOP)  AMB0421 

IMPLICIT  REAL*S(A-H,0-Z)  AMB0422 

CHARACTERS  ISTOP(l)  AMB0423 
COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G3,G9,G10,G11,G12,G13,G14,G15,AMB0424 

1             G16,G17,G18,G19,G20  AMB0425 

COMMON  /PAR/CO,ENO,EM1,D,SIGMA,TLIM,DRO,ELO,QO,TO,FACT,ALOGF,  AMB 0426 

1  DPSIO, DTMAX, DETAO, ETALIM, XSI , XSF  AMB0427 
COMMON  /NPAR/NPHI,IPAR,NP,NR,NX,NXS,NS,NSPEC,NS1,NS2,NTAU0,NETA0,  AMB0423 

1             NAMB,NCASE,ICASE,IFAN  AMB0429 

COMMON  /GE0M/APF,PAI,PAI2,W,SW,CW,BETA,SBETA,CBETA,PSI1,SPSI1,  AMB 0  43  0 

1  CPSI1,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,A0,RF,XF,YF,ZF,AMB0431 

2  PHISOF,PHIF,SPHIF,CPHIF,DYMIN,RMIN,XS,DIST,X0,Y0,20,  AMB0432 
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COMMON 
COMMON 


DYO , DEC  PSIN, ST1 , CT1 , OMEGX, OMEGY, OMEGZ, XSVC  21 ) 
/EPSIL/EPSETA, EPST, EPSR 


ETAKXTC5),PHIEXT(5), 

FEXTC5),WEXT(5), 
ETAMAX(5),PSIMAX(5), 


PHIMAX(5),WMAX(5) 


/EXTREM/TEXT(5),ETAEXT(5) 

1  PSIEXT(5),EMEXT(5) 

2  TMAX(5),ETAKMX(5) 

3  EMMAX(5),FMAX(5), 

4  RFMAX(5),PHIFMX(5) 
COMMON  /SOFPR/C, DSUMF, DSUMD, T, ETA, DETA, SUM, DSUM, SUMU, DSUMU 
COMMON  /SUMS/SUMFC5),SUMD(5),SUMU2(5) 

COMMON  /COUNTS/ICONTC, ICONTT, ICNTOT, ICNTMX, IQT0TC5) , ISHADC5) 
COMMON  /SPEC/WAV,XCC5),WC(5),WRC(5),XNAME(5),QFC(5),QDCC5), 
1  QU2C(5),FLUXCC5),0MEGA(5),FLXU2CC5),URMSCC5) 

PRINT  l,ISTOP 
1     FORMAT(///1X,2H**,2X,30A4,2X,2H**,///) 
PRINT  71 , NS, NP, NR, NX, ICONTC, ICONTT 

71  FORMATC1X, »NS, NP, NR, NX, ICONTC, ICONTT= ' ,616/) 
IF(NS.GT.NSPEC)  NS=1 

PRINT  72,RF,PHIF*DEG,PHIS0F*DEG,W*DEG,BETA*DEG 

72  FORMATC/1X, 'RF, PHIF, PHISOF, W, BETA= ' , D14 . 5, 4F10 . 3/) 
PRINT  73,C,T,TLIM,ETA 

73  FORMATC/1X,  'C, T, TLIM, ETA= ' , 4D14 . 5/) 

PRINT  74, DSUM, SUM, DSUMF, SUMF(NS) , SUMD(NS) , QDC(NS) , QFC(NS) , 
1  FLUXC(NS),OMEGA(NS) 

74  FORMATC IX, ■ DSUM, SUM, DSUMF, SUMF(NS) , SUMDC  NS)  =  ■ , 5D14 . 5/ 

1        IX, •QDC(NS),QFC(NS),FLUXC(NS),0MEGA(NS)=,,4D14.5/) 

XX=-1.D0 

YY=DSQRT(XX)+1.D0 

STOP 

END 

SUBROUTINE  FLUX 

IMPLICIT  REAL*8(A-H,0-Z) 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15 
1  G16,G17,G18,G19,G20 

COMMON  /PAR/CO, ENO, EMI, D, SIGMA, TLIM,DR0, EL 0,Q0, TO, FACT, ALOGF, 
1  DPSIO,DTMAX,DETAO,ETALIM,XSI,XSF 

COMMON  /NPAR/NPHI,IPAR,NP,NR,NX,NXS,NS,NSPEC,NS1,NS2,NTAU0,NETA0, 
1  NAMB,NCASE,ICASE,IFAN 

COMMON  /GE0M/APF,PAI,PAI2,W,SW,CW,BETA,SBETA,CBETA,PSI1,SPSI1, 

1  CPSIl,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,AO,RF,XF,YF,ZF 

2  PHISOF, PHI F,SPHIF,CPHIF,DYMIN,RMIN,XS,DIST,XO,YO,ZO, 

3  DYO , DEG, PSIN, ST1 , CT1 , OMEGX, OMEGY, OMEGZ, XSV( 21 ) 
COMMON  /EPSIL/EPSETA, EPST, EPSR 

COMMON  /EXTREM/TEXT ( 5 ) , ETAEXTC  5 ) , ETAKXT ( 5 ) , PHI EXT( 5 ) , 

1  PSIEXT(5),EMEXT(5),FEXT(5),WEXT(5), 

2  TMAXC5) , ETAKMXC5) , ETAMAXC5) , PSIMAX(5) , 

3  EMMAX(5),FMAX(5), 

4  RFMAX(5),PHIFMX(5),PHIMAX(5),WMAX(5) 
COMMON  /SOFPR/C, DSUMF, DSUMD, T, ETA, DETA, SUM, DSUM, SUMU, DSUMU 
COMMON  /COUNTS/ ICONTC, ICONTT, ICNTOT, ICNTMX, I QTOT( 5), I SHADC 5) 
COMMON  /SPEC/WAV,XC(5),WC(5),WRC(5),XNAME(5),QFC(5),QDC(5), 

1  QU2C(5),FLUXC(5),0MEGA(5),FLXU2C(5),URMSC(5) 

SUMD(5),SUMU2(5) 


HERE  IS  NOT  WRITTEN  FOR  ZO.NE.O.') 


COMMON  /SUMS/SUMFC5); 

ELO=SIGMA*RF*ENO 

IF(ZO.NE.O.) 
1CALL  SOF( 'THE  SCHEME 

YYO=(YO-AO)/XO 

PCHECK=DATAN(YYO) 

IFCPCHECK.GT.PSIF-1.D-4.0R.PCHECK.LT.-1.D-4) 
1CALL    SOFCFLUX    RECEIVING    POINT    WITHIN    FAN    OR 

SPHIF=DSIN(PHIF) 

CPHIF=DCOS(PHIF) 

XF=RF*CPSIF 

YF=APF*CPHIF 

ZF=APF*SPHIF 

TBETA=ZF/(YF-YO) 

BETA=DATAN(TBETA) 

IFCDABS(BETA) . GT 

SBETA=DSIN(BETA) 

CBETA=DCOS(BETA) 

DIST=DSQRT((XF-X0)**2+(YF-Y0)x*2+(ZF-Z0)**2) 

CW=(XF-XO)/DIST 

SW=DSQRT(l.DO-CW*X2) 


WITHIN  SPACECRAFT') 


PAI2)  CALL  SOFCBETA.GT.PAI/2') 
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W=PAI2-DATAN(CW/SW) 
OMEGX=CW 
OMEGY=SW*CBETA 
0MEGZ=SW*SBETA 
CALL  LIMIT 


NS2 


DO  20  NS=NS1 

SUMF(NS)=0. 

SUMU2(NS)=0. 

SUMD(NS)=0. 

FEXT(NS)=0. 

CALL  SUMT 

SUMF(NS)=SUM 

SUMU2(NS)=SUMU 

QFC(NS)=SUMF(NS)/FACT 

QU2C(NS)=SUMU2(NS)/FACT 

FEXT(NS)=FEXT(NS)/FACT 

CALL  FAN(TEXT(NS),PSIEXTCNS),PHIEXT(NS)) 

IF(PSIEXT(NS).LT.PSIF-1.D-10)  CALL  SOF( • PSIEXT(NS) . LT . PSIF' ) 

IFCPSIEXT(NS).GT.PSIl)    PSIEXT(NS)=PSI1 

PSIO=PSIEXT(NS) 

T=TEXT(NS) 

CALL  MATCH(T,PSIO,EM,TETA) 

EMEXT(NS)=EM 

WEXTCNS)=W 

IQTOT(NS)=IQTOTCNS)+ICONTT 
20    CONTINUE 

RETURN 

END 

SUBROUTINE  LIMIT 

IMPLICIT  REAL*8(A-H,0-Z) 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15, 
1  G16,G17,G18,G19,G20 

COMMON  /PAR/CO, ENO, EMI, D, SIGMA, TLIM,DRO, EL 0,Q0, TO, FACT, ALOGF, 
1  DPSIO,DTMAX,DETAO,ETALIM,XSI,XSF 

COMMON  /NPAR/NPHI,IPAR,NP,NR,NX,NXS,NS,NSPEC,NS1,NS2,NTAU0,NETA0, 
1  NAMB,NCASE,ICASE,IFAN 

COMMON  /GE0M/APF,PAI,PAI2,W,SW,CW,BETA,SBETA,CBETA,PSI1,SPSI1, 

1  CPSI1,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,A0,RF,XF,YF,ZF, 

2  PHISOF,PHIF,SPHIF,CPHIF,DYMIN,RMIN,XS,DIST,X0,Y0,Z0, 

3  DY0,DEG,PSIN,ST1,CT1,0MEGX,0MEGY,0MEGZ,XSV(21) 
COMMON  /EPSIL/EPSETA,EPST,EPSR 

COMMON  / EXTREM/TEXT ( 5 ) , ETAEXT ( 5 ) , ETAKXT ( 5 ) , PHI  EXT ( 5 ) , 
PSIEXT(5),EMEXT(5),FEXT(5),WEXT(5), 
TMAX(5),ETAKMX(5),ETAMAX(5),PSIMAX(5), 
EMMAX(5),FMAX(5), 
RFMAX(5),PHIFMX(5),PHIMAX(5),WMAX(5) 

COMMON  /SPEC/WAV,XC(5),WC(5),WRC(5),XNAME(5),QFC(5),QDC(5), 
QU2C(5),FLUXC(5),0MEGA(5),FLXU2C(5),URMSC(5) 

AAA=(CW/CPSI1)**2-1 .DO 

IF(AAA.LT.l.D-lO)  GO  TO  1 

TPSI1=SPSI1/CPSI1 

AP1=A0+XF*TPSI1 

BBB=2.D0X(AP1*CW*TPSI1-SW*APF*(CBETA*CPHIF+SBETA*SPHIF)) 

CCC=AP1**2-APF*X2 

DDD=BBB**2-4.D0*AAA*CCC 

TLIM=(-BBB+DSQRT(DDD))/(2.D0*AAA) 

TLIM=TLIM/RF 

RETURN 

CONTINUE 

TLIM=1.D  55 

RETURN 

END 

SUBROUTINE  SUMT 

IMPLICIT  REAL*8(A-H,0-Z) 

COMMON  /GAMA/G, Gl, G2, G3, G4, G5,G6,G7,G8,G9,G10,G11,G1 2, G13,G1 4, G15, 

G16,G17,G18,G19,G20 
COMMON  /PAR/CO, ENO, EMI, D, SIGMA, TLIM,DRO, EL 0,Q0, TO, FACT, ALOGF, 

DPSIO,DTMAX,DETAO,ETALIM,XSI,XSF 
COMMON  /NPAR/NPHI,IPAR,NP,NR,NX,NXS,NS,NSPEC,NS1,NS2,NTAU0,NETA0, 

NAMB,NCASE,ICASE,IFAN 
COMMON  /GE0M/APF,PAI,PAI2,W,SW,CW,BETA,SBETA,CBETA,PSI1,SPSI1, 
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COMMON 
COMMON 


COMMON 
COMMON 
COMMON 


CPSI1,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,A0,RF,XF,YF,ZF, 
PHISOF,PHIF,SPHIF,CPHIF,DYMIN,RMIN,XS,DIST,X0,Y0,Z0, 
DYO , DEG, PSIN, ST1 , CT1 , OMEGX, OMEGY, OMEGZ, XSV( 21) 
/EPSIL/EPSETA,EPST,EPSR 

/EXTREM/TEXT(5),ETAEXT(5),ETAKXT(5),PHIEXT(5), 
PSIEXT(5),EMEXT(5),FEXT(5),WEXT(5), 
TMAXC  5) , ETAKMXC  5) , ETAMAXC  5) , PSIMAXC5) , 
EMMAX(5),FMAX(5), 

RFMAXC5),PHIFMX(5),PHIMAX(5),WMAXC5) 
/SOFPR/CC, DSUMF, DSUMD, T, ETA, DETA, SUM, DSUM, SUMU, DSUMU 
/COUNTS/ICONTC, ICONTT, ICNTOT, ICNTMX, IQT0TC5) , ISHADC5) 
/SPEC/WAV,XC(5),WCC5),WRC(5),XNAME(5),QFC(5),QDCC5), 

QU2CC5),FLUXC(5),0MEGA(5),FLXU2C(5),URMSC(5) 
/SUMS/SUMFC5),SUMD(5),SUMU2(5) 

OF  FLUX  ARRIVING  ALONG  A  SINGLE  RAY 


10 


15 


111 


COMMON 
INTEGRATION 

DT=DPSI0 

PSIN=PSIF 

ETA1=0. 

ETA3=0. 

FETA4=0. 

FETAU4=0. 

T  =  0. 

SUM=0. 

SUMU=0. 

ICONTT=0 

ICONTT=ICONTT+l 

PSIL=PSIN 

DT2=DT/2.D0 

DT6=DT/6.D0 

T1=T+DT2 

T2=T+DT 

FETA1=FETA4 

FETAU1=FETAU4 

CALL  PATHK(T1,ETAK1) 

CAL L  FETACT1 , ETA1 , ETAK1 , GT2, FETA2, FETAU2) 

FETA3=FETA2 

FETAU3=FETAU2 

CALL  PATHK(T2,ETAK3) 

CALL  FETACT2,  ETA3,  ETAK3,GT4,  FETA4, FETAU4) 

DETA=DT*GT4 

DSUM=DT6*(FETA1+2.D0*(FETA2+FETA3)+FETA4) 

DSUMU=DT6*(FETAU1+2.D0*(FETAU2+FETAU3)+FETAU4) 

T=T+DT 

ETA=ETA3 

ETAK=ETAK3 

SUM=SUM+DSUM 

SUMU=SUMU+DSUMU 

IF(FEXTCNS) .GT.FETA4)  GO  TO  10 

FEXT(NS)=FETA4 

TEXT(NS)=T 

ETAEXT(NS)=ETA 

ETAKXT(NS)=ETAK 

CONTINUE 
STEP  CONTROL  (DT) 

CALL  FAN(T,PSI,PHI) 

IF(PSI.LT.PSIF-l.D-lO)  CALL  SOF( ' PSI . LT . PSIF' ) 

IF(PSI.GT.PSIl)  PSI=PSI1 

PSIN=PSI 

DPSI=PSIN-PSIL 

DTP=DT*(DPSI0/(DPSI+1.D-10)) 

DTE=DT*(DETA0/(DETA+1.D-10)) 

DT1=1.2D0XDT 

DT=DMIN1(DTP,DTE,DT1,DTMAX) 

IF(DT.LE.O.)  CALL  SOF( 'COMPUTED  DT  NEGATIVE*) 

CONTINUE 

IFCIPAR.LT. 1) 
1PRINT  1 11,NR,NP,T,PSI*DEG,  PHI *DEG,  ETA,  ETAK,  SUM,  DSUM/ (SUM+1 

FORMATdX,  lNR,NP,T,PSI,PHI  =  ,,2I3,3D12.3/ 
1        IX,  'ETA,ETAK,SUM,ERRR=',4D12.3) 

IFCICONTT.GT. ICNTMX) 
1CALL  SOF( 'ICONTT  TOO  LARGE') 

IF(IC0NTT.LE.2)  GO  TO  1 


D-20) 
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IFCETA+ETAK.GT.ETALIM)  GO  TO  100  AMB0649 

IFCT.GT.50.D0  .OR.  T*RF.GT.A0)  GO  TO  100  AMB0650 

IFCSUM.EQ.O.)  GO  TO  1  AMB0651 

ERR=(DSUM/SUM)/DT  AMB0652 

IF(ERR.GT.EPST)  GO  TO  1  AMB0653 

100   CONTINUE  AMB0654 

SUM=SUM*ELO  AMB0655 

SUMU=SUMU*ELO  AMB0656 

RETURN  AMB0657 

END  AMB0658 

SUBROUTINE  FETACT, ETAIK, ETAK,GT, FET, FETU2)  AMB0659 

IMPLICIT  REAL*8(A-H,0-Z)  AMB0660 

REALX8  MU1,MU2  AMB0661 
COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15,AMB0662 

1             G16,G17,G18,G19,G20  AMB0663 

COMMON  /PAR/CO,  ENO,  EMI,  D,  SIGMA,  TLIM,DR0,  EL  0,Q0,  TO,  FACT,  ALOGF,  AMB066  4 

1  DPSIO,DTMAX,DETAO,ETALIM,XSI,XSF  AMB0665 
COMMON  /NPAR/NPHI,IPAR,NP,NR,NX,NXS,NS,NSPEC,NS1,NS2,NTAU0,NETA0,  AMB0666 

1  NAMB,NCASE,ICASE,IFAN  AMB0667 
COMMON  /GE0M/APF,PAI,PAI2,W,SW,CW,BETA,SBETA,CBETA,PSI1,SPSI1,     AMB0668 

1  CPSI1,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,A0,RF,XF,YF,ZF,AMB066  9 

2  PHISOF,PHIF,SPHIF,CPHIF,DYMIN,RMIN,XS,DIST,X0,Y0,Z0,  AMB0670 

3  DYO,DEG,PSIN,ST1,CT1,OMEGX,OMEGY,OMEGZ,XSV(21)  AMB0671 
COMMON  /EPSIL/EPSETA,EPST,EPSR  AMB0672 
COMMON  /EXTREM/TEXT(5),ETAEXT(5),ETAKXT(5),PHIEXT(5),  AMB0673 

PSIEXTC5),EMEXT(5),FEXT(5),WEXT(5),  AMB067  4 

TMAX(5),ETAKMX(5),ETAMAX(5),PSIMAX(5),  AMB0675 

EMMAX(5),FMAX(5),  AMB0676 

RFMAX(5),PHIFMX(5),PHIMAX(5),WMAX(5)  AMB0677 

COMMON  /SPEC/WAV,XC(5),WC(5),WRC(5),XNAME(5),QFC(5),QDC(5),  AMB0678 

1             QU2C(5),FLUXC(5),0MEGA(5),FLXU2C(5),URMSC(5)  AMB0679 

COMMON  /AMBIEN/ENA,UA,PSIA,PHIA,HA(3),WA(3),  AMB0680 

1               UAX,UAY,UAZ,AA,BA,CA,RA,XA,YA,ZA, SHADOW  AMB0681 

LOGICAL  SHADOW  AMB0682 

COMMON  /NAGESH/PIK,UIK,UIKX,UIKY,UIKZ  AMB0683 

ETAIK  =  0.  AMB0684 

IF(SHADOW)  GO  TO  1  AMB0685 

K=l  AMB0686 

I=NS  AMB0687 

CALL  FAN(T,PSI,PHI)  AMB0688 

IFCPSI.LT.PSIF-l.D-10)  CALL  SOF( ' PSI . LT . PSIF' )  AMB0689 

IF(PSI.GT.PSIl)  PSI=PSI1  AMB0690 

PSI0=PSI  AMB0691 

CALL  MATCH(T,PSIO,EM,TETA)  AMB0692 

SPSI=DSIN(PSI)  AMB0693 

CPSI=DCOS(PSI)  AMB0694 

SPHI=DSIN(PHI)  AMB0695 

CPHI=DCOS(PHI)  AMB0696 

ST=DSIN(TETA)  AMB0697 

CT=DCOS(TETA)  AMB0698 

G0REM=1.D0+G1*EM**2  AMB0699 

TERMN=G0REM**G6  AMB0700 

U=EM*C0/DSQRT(GOREM)  AMB0701 

UX=U*CT  AMB0702 

UY=U*ST*CPHI  AMB0703 

UZ=U*ST*SPHI  AMB0704 

COLLISION  AMB0705 

MU1=WC(I)/(WC(I)+WA(K))  AMB0706 

MU2=1.D0-MU1  AMB0707 

UMX=MU1*UX+MU2*UAX  AMB0708 

UMY=MU1*UY+MU2*UAY  AMB0709 

UMZ=MU1*UZ+MU2*UAZ  AMB0710 

DOTUM=OMEGX*UMX+OMEGY*UMY+OMEGZ*UMZ  AMB0711 

URX=UX-UAX  AMB0712 

URY=UY-UAY  AMB0713 

URZ=UZ-UAZ  AMB0714 

UR=DSQRT(URX*x2+URY**2+URZ**2)  AMB0715 

DET=D0TUM**2+(MU2*UR)**2-(UMX*X2+UMY**2+UMZ**2)  AMB  07 16 

IFCDET.LT.O.)  GO  TO  1  AMB0717 

DET1=DSQRT(DET)  AMB0718 

UIK1=-D0TUM+DET1  AMB0719 

UIK2=-D0TUM-DET1  AMB0720 
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IFCUIK2.GT.0.)  CALL  SOFCDOUBLE  COLLISION  OPTION  NOT  PROGRAMMED    AMB0721 

1  YET')  AMB0722 

UIK=UIK1  AMB0723 

IFCUIK.LE.O.)  GO  TO  1  AMB0724 

UIKX=-OMEGX*UIK  AMB0725 

UIKY=-OMEGYxUIK  AMB0726 

UIKZ=-OMEGZXUIK  AMB0727 

CDEL=(D0TUM+UIK)/(MU2XUR)  AMB0728 

IFCCDEL.LE.O.)  CALL  SOFCCDEL  NEGATIVE  NOT  PROGRAMMED  YET')        AMB0729 

IF(CDEL-l.D-lO.GT.l.DO)  AMB0730 

1CALL  SOFCCDEL  CCOSCDELTA))  CANNOT  BE  GT  .  1 .  »  )  AMB0731 

PIK=(UIK/(MU2*UR))*X2/(4.D0*PAI*CDEL)  AMB0732 

IF  (PIK.LT.O.)  CALL  SOFC ■ PIK . LT . 0 ' )  AMB0733 

FET=(UR/UA)*PIK/TERMN  AMB0734 

UREL=DSQRT((UX-UIKX)XX2+(UY-UIKY)XX2+(UZ-UIKZ)XX2)  AMB0735 

GT  =  ELOX(UREL/UIKVTERMN  AMB0736 

CALL  PATHIKCT,ETAIK)  AMB0737 

POWER=ETAIK+ETAK-ALOGF  AMB0738 

EFACT=0.  AMB0739 

IFC POWER. LT. 60. DO )EFACT=DEXP( -POWER)  AMB0740 

FET=FETXEFACT  AMB0741 

FETU2  =  FETXUIKXX2  AMB0742 

IFCEM.LT.O.)  CALL  SOFC ' EM . LT . 0 • )  AMB0743 

FETU2=FETXEM  AMB0744 

RETURN  AMB0745 

CONTINUE  AMB0746 

FET=0.  AMB0747 

FETU2=0.  AMB0748 

GT=0.  AMB0749 

RETURN  AMB0750 

END  AMB0751 

SUBROUTINE  PATHIK(TC, ETAIK)  AMB0752 

IMPLICIT  REALX8(A-H,0-Z)  AMB0753 

REALX8  MU1,MU2  AMB0754 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15,AMB0755 

1              G16,G17,G18,G19,G20  AMB0756 

COMMON  /PAR/CO, ENO, EMI, D, SIGMA, TLIM,DRO, EL 0,Q0, TO, FACT, ALOGF,  AMB0757 

1  DPSIO,DTMAX,DETAO,ETALIM,XSI,XSF  AMB0758 
COMMON  /NPAR/NPHI,IPAR,NP,NR,NX,NXS,NS,NSPEC,NS1,NS2,NTAU0,NETA0,  AMB0759 

1  NAMB,NCASE,ICASE,IFAN  AMB0760 
COMMON  /GE0M/APF,PAI,PAI2,W,SW,CW,BETA,SBETA,CBETA,PSI1,SPSI1,     AMB0761 

1  CPSI1,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,A0,RF,XF,YF,ZF,AMB0762 

2  PHISOF,PHIF,SPHIF,CPHIF,DYMIN,RMIN,XS,DIST,X0,Y0,Z0,  AMB0763 

3  DY0,DEG,PSIN,STl,CTl,OMEGX,OMEGY,OMEGZ,XSV(21)  AMB0764 
COMMON  /EPSIL/EPSETA,EPST,EPSR  AMB0765 
COMMON  /EXTREM/TEXT(5),ETAEXTC5),ETAKXT(5),PHIEXT(5),  AMB0766 

1  PSIEXT(5),EMEXT(5),FEXT(5),WEXT(5),  AMB0767 

2  TMAX(5),ETAKMX(5),ETAMAX(5),PSIMAX(5),  AMB0768 

3  EMMAX(5),FMAX(5),  AMB0769 

4  RFMAX(5),PHIFMX(5),PHIMAX(5),WMAX(5)  AMB0770 
COMMON  /SPEC/WAV,XC(5),WC(5),WRC(5),XNAME(5),QFC(5),QDC(5),  AMB0771 

1             QU2CC5),FLUXC(5),0MEGA(5),FLXU2C(5),URMSC(5)  AMB0772 

COMMON  /AMBIEN/ENA,UA,PSIA,PHIA,HA(3),WA(3),  AMB0773 

1               UAX,UAY,UAZ,AA,BA,CA,RA,XA,YA,ZA, SHADOW  AMB0774 

LOGICAL  SHADOW  AMB0775 

NETA=NETA0  AMB0776 

DT=TC/DBLE(NETA)  AMB0777 

DT2=DT/2.D0  AMB0778 

DT6=DT/6.D0  AMB0779 

GT4=0.  AMB0780 

T=0.  AMB0781 

ETA=0.  AMB0782 

IT=0  AMB0783 

IT=IT+1  AMB0784 

T1=T+DT2  AMB0785 

T2=T+DT  AMB0786 

GT1=GT4  AMB0787 

CALL  FT(T1,GT2)  AMB0788 

GT3=GT2  AMB0789 

CALL  FT(T2,GT<+)  AMB0790 

DETA=DT6X(GT1+2.D0*(GT2+GT3)+GT4)  AMB07  91 

T=T+DT  AMB0792 
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ETA=ETA+DETA 
IF(IT.LT.NETA)  GO  TO  1 
ETAIK=ETA 
RETURN 
END 

SUBROUTINE  FTCT,GT) 
IMPLICIT  REAL*8(A-H,0-Z) 
REAL*8  MU1,MU2 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15, 
I  G16,G17,G18,G19,G20 

COMMON  /PAR/CO, ENO, EMI, D, SIGMA, TLIM,DRO, EL 0,Q0, TO, FACT, ALOGF, 
L  DPSIO,DTMAX,DETAO,ETALIM,XSI,XSF 

COMMON  /NPAR/NPHI,IPAR,NP,NR,NX,NXS,NS,NSPEC,NSl,NS2,NTAUO,NETAO, 
L  NAMB,NCASE,ICASE,IFAN 

COMMON  /GE0M/APF,PAI,PAI2,W,SW,CW,BETA,SBETA,CBETA,PSI1,SPSI1, 

CPSIl,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,AO,RF,XF,YF,ZF, 
I  PHISOF,PHIF,SPHIF,CPHIF,DYMIN,RMIN,XS,DIST,XO,YO,ZO, 

5  DYO , DEG, PSIN, ST1 , CT1 , OMEGX, OMEGY, OMEGZ, XSV( 21 ) 

COMMON  /EPSIL/EPSETA,EPST,EPSR 
COMMON  /EXTREM/TEXTC5) , ETAEXT(5) , ETAKXT(5) , PHIEXTC5) , 

PSIEXT(5),EMEXT(5),FEXT(5),WEXT(5), 

TMAX( 5) , ETAKMX( 5) , ETAMAXC5) , PSIMAXC5) , 

EMMAX(5),FMAX(5), 

RFMAX(5),PHIFMX(5),PHIMAX(5),WMAX(5) 
COMMON  /SPEC/WAV,XC(5),WC(5),WRC(5),XNAME(5),QFC(5),QDC(5), 

QU2CC5),FLUXC(5),0MEGAC5),FLXU2C(5),URMSC(5) 
COMMON  /AMBIEN/ENA,UA,PSIA,PHIA,HA(3),WA(3), 

UAX,UAY,UAZ,AA,BA,CA,RA,XA,YA,ZA, SHADOW 
LOGICAL  SHADOW 

COMMON  /NAGESH/PIK,UIK,UIKX,UIKY,UIKZ 
K  =  l 
I  =  NS 

CALL  FAN(T,PSI,PHI) 

IF(PSI.LT.PSIF-l.D-lO)  CALL  SOF( 'PSI . LT . PSIF' ) 
IF(PSI.GT.PSIl)  PSI=PSI1 
PSI0=PSI 

CALL  MATCH(T,PSIO,EM,TETA) 
SPSI=DSIN(PSI) 
CPSI=DCOS(PSI) 
SPHI=DSIN(PHI) 
CPHI=DCOS(PHI) 
ST=DSIN(TETA) 
CT=DCOS(TETA) 
G0REM=1.D0+G1*EM**2 
TERMN=GOREM**G6 
U=EM*CO/DSQRT(GOREM) 
UX=U*CT 
UY=U*ST*CPHI 
UZ=UXST^SPHI 

UREL  =  DSQRT((UX-UIKX)**2+(UY-UIKY)**2+(UZ-UIKZ:>**2) 
GT=ELO*(UREL/UIK)/TERMN 
RETURN 
END 

SUBROUTINE  PATHKCT, ETAK) 
IMPLICIT  REAL*8(A-H,0-Z) 
COMMON  /GAMA/G, Gl, G2, G3, G4, G5,G6,G7,G8,G9,G10,G11,G12,G1 3, Gl 4,  G15, 

G16,G17,G18,G19,G20 
COMMON  /PAR/CO, ENO, EMI, D, SIGMA, TLIM,DRO, EL 0,Q0, TO, FACT, ALOGF, 

DPSIO,DTMAX,DETAO,ETALIM,XSI,XSF 
COMMON  /NPAR/NPHI,IPAR,NP,NR,NX,NXS,NS,NSPEC,NS1,NS2,NTAU0,NETA0, 

NAMB,NCASE, ICASE, IFAN 
COMMON  /GE0M/APF,PAI,PAI2,W,SW,CW,BETA,SBETA,CBETA,PSI1,SPSI1, 

CPSI1,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,A0,RF,XF,YF,ZF, 
!  PHISOF,PHIF,SPHIF,CPHIF,DYMIN,RMIN,XS,DIST,X0,Y0,Z0, 

;  DYO, DEG, PSIN, ST1,CT1, OMEGX, OMEGY, OMEGZ, XSV( 21) 

COMMON  /EPSIL/EPSETA,EPST,EPSR 
COMMON  /EXTREM/TEXT(5),ETAEXT(5),ETAKXT(5),PHIEXT(5), 

PSIEXT(5),EMEXT(5),FEXT(5),WEXT(5), 

TMAXC5),ETAKMX(5),ETAMAX(5),PSIMAX(5), 

EMMAX(5),FMAX(5), 

RFMAX(5),PHIFMX(5),PHIMAX(5),WMAX(5) 
COMMON  /COUNTS/ ICONTC, ICONTT, ICNTOT, ICNTMX, IQT0TC5) , ISHADC5) 


AMB0793 
AMB0794 
AMB0795 
AMB0796 
AMB0797 
AMB0798 
AMB0799 
AMB0800 
AMB0801 
AMB0802 
AMB0803 
AMB0804 
AMB0805 
AMB0806 
AMB0807 
AMB0808 
AMB0809 
AMB0810 
AMB0811 
AMB0812 
AMB0813 
AMB0814 
AMB0815 
AMB0816 
AMB0817 
AMB0818 
AMB0819 
AMB0820 
AMB0821 
AMB0822 
AMB0823 
AMB0824 
AMB0825 
AMB0826 
AMB0827 
AMB0828 
AMB0829 
AMB0830 
AMB0831 
AMB0832 
AMB0833 
AMB0834 
AMB0835 
AMB0836 
AMB0837 
AMB0838 
AMB0839 
AMB0840 
AMB0841 
AMB0842 
AMB0843 
AMB0844 
AMB0845 
AMB0846 
AMB0847 
AMB0848 
AMB0849 
AMB0850 
AMB0851 
AMB0852 
AMB0853 
AMB0854 
AMB0855 
AMB0856 
AMB0857 
AMB0858 
AMB0859 
AMB0860 
AMB0861 
AMB0862 
AMB0863 
AMB0864 


45 


SCRIPT   Al 

COMMON  /AMBIEN/ENA,UA,PSIA,PHIA,HA(3),WA(3),  AMB0865 

1               UAX,UAY,UAZ,AA,BA,CA,RA,XA,YA,ZA, SHADOW  AMB0866 

LOGICAL  SHADOW  AMB0867 

ETAK=0.  AMB0868 

:   DETERMINE  POINT  OF  ENTRY  OF  AMBIENT  TRAJECTORY  TO  FAN  AMB0869 

TRF=T*RF  AMB0870 

XC=XF+TRF*OMEGX  AMB0871 

YC=YF+TRF*OMEGY  AMB0872 

ZC=ZF+TRF*OMEGZ  AMB0873 

:   CHECK  SHADOW  AMB0874 

SHADOW=. FALSE.  AMB0875 

EVER=BA**2+CA**2  AMB0876 

DETS=EVER*A0XX2-(BA*ZC-CA*YC)**2  AMB0877 

IF(DETS.LE.O.)  GO  TO  2  AMB0878 

DETS1=DSQRT(DETS)  AMB0879 

TAU1=(-(BAXYC+CAXZC)+DETS1)/EVER  AMB088  0 

IFCTAU1.GT.0.)  SHADOW=.TRUE.  AMB0881 

2    CONTINUE  AMB0882 

IF(SHADOW)  GO  TO  10  AMB0883 

EVER1=A0+XCXTPSIF  AMB0884 

EVER2=BAXX2+CAXX2-(AAXTPSIF)XX2  AMB0885 

EVER3=BAXYC+CAXZC-AAXEVER1XTPSIF  AMB0886 

DET=EVER3XX2-EVER2X(YCXX2+ZCXX2-EVER1XX2)  AMB0887 

IFCDET.LE.O.)  AMB0888 

1CALL  SOFCNO  INTERSECTION  OF  AMB .  TRAJ .  WITH  LIMITING  CONE')  AMB0889 

DET1=DSQRT(DET)  AMB0890 

TAUP=(-EVER3+DET1)/EVER2  AMB0891 

TAUM=(-EVER3-DET1)/EVER2  AMB0892 

IFCTAUP.GT.O.  .AND.  TAUM.GT.O.)  AMB0893 

1CALL  SOFCTWO  POSITIVE  INTERSECTIONS  WITH  LIMITING  CONE. NOT  PERMITAMB0894 

1  IN  THIS  VERSION')  AMB0895 

TAUF=DMAX1(TAUP,TAUM)  AMB0896 

IFCTAUF.LE.O.)  AMB0897 

1CALL  SOFCNO  POSITIVE  INTERSECTION  WITH  LIMITING  CONE')  AMB0898 

XA=XC+TAUFXAA  AMB0899 

YA=YC+TAUFXBA  AMB0900 

ZA=ZC+TAUFXCA  AMB0901 

RA=DSQRT(XAXX2+(DSQRT(YAXX2+ZAXX2)-A0)XX2)  AMB0902 

TAUF=TAUF/RA  AMB0903 

NTAU  =  NTAUO  AMB0904 

DTAU=TAUF/DBLE(NTAU)  AMB0905 

ETAK=0.  AMB0906 

TAU=0.  AMB0907 

DTAU2=DTAU/2.D0  AMB0908 

DTAU6=DTAU/6 .DO  AMB0909 

GTAU4=0.  AMB0910 

ITAU=0  AMB0911 

1     ITAU=ITAU+1  AMB0912 

TAU1=TAU+DTAU2  AMB0913 

TAU2=TAU+DTAU  AMB0914 

GTAU1=GTAU4  AMB0915 

CALL  FTAU(TAU1,GTAU2)  AMB0916 

GTAU3=GTAU2  AMB0917 

CALL  FTAU(TAU2,GTAU4)  AMB0918 

DETAK=DTAU6X(GTAU1+2.D0X(GTAU2+GTAU3)+GTAU4)  AMB0919 

TAU=TAU+DTAU  AMB0920 

ETAK=ETAK+DETAK  AMB0921 

IF(ITAU.LT.NTAU)  GO  TO  1  AMB0922 

ETAK=ETAKX(SIGMAXENOXRA)  AMB0923 

RETURN  AMB0924 

10    CONTINUE  AMB0925 

ISHAD(NS)=ISHAD(NS)+1  AMB0926 

RETURN  AMB0927 

END  AMB0928 

SUBROUTINE  FTAUCTAU, GTAU)  AMB0929 

IMPLICIT  REALX8(A-H,0-Z)  AMB0930 
COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15,AMB0931 

1             G16,G17,G18,G19,G20  AMB0932 

COMMON  /PAR/CO, ENO, EMI, D, SIGMA, TLIM,DRO, EL 0,Q0, TO, FACT, ALOGF,  AMB0933 

1  DPSIO,DTMAX,DETAO,ETALIM,XSI,XSF  AMB0934 
COMMON  /NPAR/NPHI,IPAR,NP,NR,NX,NXS,NS,NSPEC,NS1,NS2,NTAU0,NETA0,  AMB0935 

1             NAMB,NCASE,ICASE,IFAN  AMB0936 
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COMMON  /GE0M/APF,PAI,PAI2,W,SW,CW,BETA,SBETA,CBETA,PSI1,SPSI1,  AMB0937 

1  CPSI1,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,A0,RF,XF,YF,ZF,AMB0938 

2  PHISOF,PHIF,SPHIF,CPHIF,DYMIN,RMIN,XS,DIST,X0,Y0,Z0,  AMB0939 

3  DY0,DEG,PSIN,ST1,CT1,0MEGX,0MEGY,0MEGZ,XSV(21)  AMB0940 
COMMON  /EPSIL/EPSETA,EPST,EPSR  AMB0941 
COMMON  /EXTREM/TEXT(5),ETAEXTC5),ETAKXT(5),PHIEXTC5),  AMB0942 

1  PSIEXT(5),EMEXT(5),FEXT(5),WEXT<5),  AMB0943 

2  TMAX(5),ETAKMX(5),ETAMAX(5),PSIMAX(5),  AMB0944 

3  EMMAXC5),FMAX(5),  AMB0945 

4  RFMAX(5),PHIFMX(5),PHIMAX(5),WMAXC5)  AMB0946 
COMMON  /SPEC/WAV, XC(5) , WCC5) , WRC( 5), XNAME( 5), QFC(5) , QDCC5) ,  AMB0947 

1             QU2CC5),FLUXC(5),0MEGA(5),FLXU2CC5),URMSC(5)  AMB0948 

COMMON  /AMBIEN/ENA,UA,PSIA,PHIA,HA(3),WA(3),  AMB0949 

1               UAX,UAY,UAZ,AA,BA,CA,RA,XA,YA,ZA, SHADOW  AMB0950 

LOGICAL  SHADOW  AMB0951 

CALL  FANT(TAU,PSI,PHI)  AMB0952 

IF(PSI.LT.PSIF-l.D-lO)  CALL  SOFC • PSI . LT . PSIF' )  AMB0953 

IF(PSI.GT.PSIl)  PSI  =  PSI1  AMB0954 

PSI0=PSI  AMB0955 

CALL  MATCH(T,PSIO,EM,TETA)  AMB0956 

SPSI=DSIN(PSI)  AMB0957 

CPSI=DCOS(PSI)  AMB0958 

SPHI=DSIN(PHI)  AMB0959 

CPHI=DCOS(PHI)  AMB0960 

ST=DSIN(TETA)  AMB0961 

CT=DCOS(TETA)  AMB0962 

GOREM=1.DO+G1*EM**2  AMB0963 

TERMN=G0REM**G6  AMB0964 

U=EM*CO/DSQRT(GOREM)  AMB0965 

UREL=DSQRT((CT*U-UAX)*X2+(ST*CPHI*U-UAY)**2+(STXSPHI*U-UAZ)*X2)  AMB0966 

GTAU=UREL/(UA*TERMN)  AMB0967 

RETURN  AMB0968 

END  AMB0969 

SUBROUTINE  FANCT, PSI , PHI )  AMB0970 

IMPLICIT  REALX8(A-H,0-Z)  AMB0971 

COMMON  /PAR/CO, ENO, EMI, D, SIGMA, TLIM,DR0, EL 0,QO, TO, FACT, ALOGF,  AMB0972 

1            DPSIO,DTMAX,DETAO,ETALIM,XSI,XSF  AMB0973 

COMMON  /GE0M/APF,PAI,PAI2,W,SW,CW, BETA, SBETA,CBETA, PSI 1,SPSI1,  AMB0974 

1  CPSI1,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,A0,RF,XF,YF,ZF,AMB097  5 

2  PHISOF,PHIF,SPHIF,CPHIF,DYMIN,RMIN,XS,DIST,XO,YO,ZO,  AMB0976 

3  DYO,DEG,PSIN,ST1,CT1,OMEGX,OMEGY,OMEGZ,XSV(21)  AMB0977 
COMMON  /POINT/XP,YP,XCOR,YCOR  AMB0978 

C   RING  FAN  GEOMETRY.    FAN  CORNER  IS  AT  ( 0 , A0*COS( PHI ) , AO*SIN( PHI) ) .  AMB0979 

C   RF   --   RADIAL  DISTANCE  ON  LIMITING  CHARACTERISTIC  OF  POINT  OF  AMB0980 

C          ENTRANCE  OF  RAY.  AMB0981 

C   DIRECTION  COSINES  OF  RAY:  OMEGX, OMEGY, OMEGZ  AMB0982 

TRF=T*RF  AMB0983 

X=XF+TRF*OMEGX  AMB0984 

Y=YF+TRF*OMEGY  AMB0985 

Z=ZF+TRF*OMEGZ  AMB0986 

DY=DSQRT(Y*Y+Z*Z)-AO  AMB0987 

IFCDABS(DY).LE.l.D-lOXAO)  DY=1.D-10*A0  AMB0988 

IFCDY.LT.O.)  AMB0989 

1CALL  SOFC'POINT  X,Y,X  CANNOT  BE  CLOSER  TO  X-AXIS  THAN  RADIUS  AO')  AMB0990 

YY=X/DY  AMB0991 

PSI=PAI2-DATAN(YY)  AMB0992 

PHI=DATAN(Z/Y)  AMB0993 

XP  =  XCOR+X  AMB0994 

YP=A0+DY  AMB0995 

RETURN  AMB0996 

END  AMB0997 

SUBROUTINE  FANTCTAU , PSI , PHI)  AMB0998 

IMPLICIT  REAL*8(A-H,0-Z)  AMB0999 

COMMON  / PAR/CO, ENO, EMI, D, SI GMA, TL IM, DRO , EL 0,Q0, TO, FACT, ALOGF,  AMB100O 

1            DPSIO,DTMAX,DETAO,ETALIM,XSI,XSF  AMB1001 

COMMON  /GE0M/APF,PAI,PAI2,W,SW,CW,BETA,SBETA,CBETA,PSI1,SPSI1,  AMB10  02 

1  CPSI1,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,A0,RF,XF,YF,ZF,AMB100  3 

2  PHISOF,PHIF,SPHIF,CPHIF,DYMIN,RMIN,XS,DIST,X0,Y0,Z0,  AMB1004 

3  DY0,DEG,PSIN,ST1,CT1, OMEGX, OMEGY, OMEGZ, XSVC21)  AMB1005 
COMMON  /AMBIEN/ENA,UA,PSIA,PHIA,HA(3),WA(3),  AMB1006 

1               UAX,UAY,UAZ,AA,BA,CA,RA,XA,YA,ZA, SHADOW  AMB1007 

COMMON  /POINT/XP,YP,XCOR,YCOR  AMB1008 


47 


)       SCRIPT    Al 

LOGICAL  SHADOW  AMB1009 

C   RING  FAN  GEOMETRY.    FAN  CORNER  IS  AT  ( 0, A0*COS(PHI) , AO*SIN(PHI) ) .    AMB1010 

C   RA   ~   RADIAL  DISTANCE  ON  LIMITING  CHARACTERISTIC  OF  POINT  OF  AMB1011 

C           ENTRANCE  OF  RAY.  AMB1012 

C   DIRECTION  COSINES  OF  RAY:  -AA,-BA,-CA  AMB1013 

TRA=TAU*RA  AMB1014 

X=XA-TRA*AA  AMB1015 

Y=YA-TRA*BA  AMB1016 

Z=ZA-TRA*CA  AMB1017 

DY=DSQRT(Y*Y+Z*Z)-AO  AMB1018 

IF(DABSCDY).LE.1.D-10*A0)  DY=1.D-10*A0  AMB1019 

IFCDY.LT.O.)  AMB1020 

1CALL  SOFCPOINT  X,Y,X  CANNOT  BE  CLOSER  TO  X-AXIS  THAN  RADIUS  A0«)  AMB1021 

YY=X/DY  AMB1022 

PSI=PAI2-DATAN(YY)  AMB1023 

PHI=DATAN(Z/Y)  AMB1024 

XP=XCOR+X  AMB1025 

YP=AO+DY  AMB1026 

RETURN  AMB1027 

END  AMB1028 

SUBROUTINE  HMSET  AMB1029 

C   SUBROUTINE  NUMBER   20  AMB1030 

IMPLICIT  REAL*8(A-H,0-Z,$)  AMB1031 

REAL*8  KAPAOB,MHINV,MINVO,M,MF,M1,M2,M3,NORM,MEXIT,LAMDOB  AMB1032 

COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15,AMB1033 

1              G16,G17,G18,G19,G20  AMB1034 

COMMON  /PAR/CO, ENO, EMI, D, SIGMA, TLIM, DRO, ELO, QO, TO, FACT, ALOGF,  AMB1035 

1             DPSIO,DTMAX,DETAO,ETALIM,XSI,XSF  AMB1036 

COMMON  /GE0M/APF,PAI,PAI2,W,SW,CW,BETA,SBETA,CBETA,PSI1,SPSI1,     AMB1037 

1  CPSI1,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,A0,RF,XF,YF,ZF,AMB1038 

2  PHISOF,PHIF,SPHIF,CPHIF,DYMIN,RMIN,XS,DIST,X0,Y0,Z0,  AMB1039 

3  DYO,DEG,PSIN,ST1,CT1,OMEGX,OMEGY,OMEGZ,XSV(21)  AMB1Q40 
COMMON  /GRP/DMINV,MHINV(101),HMV(101)  AMB1041 
COMMON  /IGRP/KHM  AMB1042 

C   A  ROUTINE  FOR  THE  C+  DERIVATIVE  DUE  TO  RING  SYMMETRY  CGRP).  AMB1043 

MEXIT  =  EM1  AMB1044 

KHM=51  AMB1045 

IF(KHM.GT.lOl)  CALL  SOF('2001')  AMB1046 

MINV0=1.D0/MEXIT  AMB1047 

DMINV=MINVO/DBLE(KHM-]J  AMB1048 

M=MEXIT  AMB1049 

SUM=0.  AMB1050 

KHM1=KHM-1  AMB1051 

DO  1  I=1,KHM1  AMB1052 

MF=M  AMB1053 

MHINV(I)=MINV0-DBLE(I-1)*DMINV  AMB1054 

M=1.D0/MHINV(I)  AMB1055 

DM=M-MF  AMB1056 

M1=M-DM  AMB1057 

M2=M-DM/2.D0  AMB1058 

M3=M  AMB1059 

CALL  MFUNC(M1,F1,ETALF1,TETA1)  AMB1060 

CALL  MFUNC(M2,F2,ETALF2,TETA2)  AMB1061 

CALL  MFUNCCM3,F3,ETALF3,TETA3)  AMB1062 

SUM  =  SUM+DM*(Fl  +  <+.D0*F2+F3)/6  .DO  AMB1063 

ETALF=ETALF3  AMB1064 

TETA=TETA3  AMB1065 

PSI=TETA+DASIN(1.D0/M)  AMB1066 

NORM=((3.D0-G)/4.DO)X(MXX2-l.D0)**O.75D0/  AMB1067 

1      (DSIN(PSI)*(1.D0+G1*M**2)**G14)  AMB1068 

HM=SUM*NORM  AMB1069 

HMV(I)=HM  AMB1070 

G0REM=1 .D0+G1*M**2  AMB1071 

G0R=M**2-1 -DO  AMB1072 

DELT0B=0.5D0*DSQRT(GOR)*(l.D0/(MEXIT*ETALF)  AMB107  3 

1       +DSIN(TETA)/M)/DSIN(PSI)+G15XHM/2.D0  AMB1074 

EPSI0B=DELT0B/DSQRT(GOR)-DSIN(TETA)/(M*DSIN(PSI))  AMB107  5 

KAPA0B=1.D0  AMB1076 

IF(DABS(PAI2-TETA) .GT.l.D-6)  AMB1077 

1KAPA0B=DTAN(TETA)*EPSI0B  AMB1078 

LAMD0B=EPSI0B-DELT0B*GOREM/(GOR*DSQRT(GOR))  AMB107  9 

PRINT  11,I,M,HM,TETA*DEG,PSI*DEG  AMB1080 
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AMB  SCRIPT        Al 

11  F0RMATC/1X,'  I,M,HM,TETA,PSI=M5,5D12.4)  AMB1081 

PRINT  12,DELTOB,EPSIOBXDEG,KAPAOB*DEG,LAMDOB*DEG  AMB1082 

12  FORMATC  IX, ' DELTOB, EPSIOB, KAPAOB, LAMD0B= ', 5X, 5D12 . 4)  AMB1083 

I  CONTINUE  AMB1084 
MHINV(KHM)=0.  AMB1085 
HMV(KHM)=l.DO  AMB1086 
RETURN  AMB1087 
END  AMB1088 
SUBROUTINE  MFUNCCM, F, ETALF, TETA)  AMB1089 

C   SUBROUTINE  NUMBER   21  AMB1090 

IMPLICIT  REALX8(A-H,0-Z,$)  AMB1091 

REALX8  NU,NUFUNC,M,MEXIT,MD,MDD  AMB1092 
COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15,AMB1093 

1  G16,G17,G18,G19,G20  AMB1094 

COMMON  /PAR/CO, ENO, EMI, D, SIGMA, TLIM,DRO, EL 0,QO, TO, FACT, ALOGF,  AMB1095 

1  DPSIO,DTMAX,DETAO,ETALIM,XSI,XSF  AMB1096 

COMMON  /GE0M/APF,PAI,PAI2,W,SW,CW,BETA,SBETA,CBETA,PSI1,SPSI1,     AMB1097 

1  CPSI1,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,A0,RF,XF,YF,ZF,AMB1098 

2  PHISOF,PHIF,SPHIF,CPHIF,DYMIN,RMIN,XS,DIST,XO,YO,ZO,  AMB1099 

3  DY0,DEG,PSIN,STl,CTl,OMEGX,OMEGY,OMEGZ,XSV(21)  AMBllOO 
C  AMB1101 

QF(MDD)=l.DO/DSQRT(MDDXX2-l.DO)  AMB1102 

NUFUNC(MD)=-G5XDATAN(G5XQF(MD))+DATANCQFCMD))  AMB1103 

C  AMB1104 

MEXIT=EM1  AMB1105 

NU=NUFUNC(M)  AMB1106 

TETA=NUFUNC(MEXIT)+PAI2-NU  AMB1107 

GOREM=1.DO+G1XMXX2  AMB1108 

GOR=M*X2-1.DO  AMB1109 

F=(MXX2)x(G0REMxxG13)XDSIN(TETA)/GORxxl.25D0  AMB1110 

GOREM1=1.DO+G1XMEXITXX2  AMB1111 

G0R1=MEXITXX2-1.D0  AMB1112 

ETALF=((GOREM/GOREMDXXGl<i)X((GORl/GOR)xx0.25DO)  AMB1113 

RETURN  AMB1114 

END  AMB1115 

SUBROUTINE  HINTER(M,H)  AMB1116 

C   SUBROUTINE  NUMBER   22  AMB1117 

IMPLICIT  REAL*8(A-H,0-Z,$)  AMB1118 

REALX8  MINV,M,MEXIT,MHINV  AMB1119 
COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15,AMB1120 

1  G16,G17,G18,G19,G20  AMB1121 

COMMON  /PAR/CO, ENO, EMI, D, SIGMA, TLIM,DRO, EL 0,QO, TO, FACT, ALOGF,  AMB1122 

1  DPSIO,DTMAX,DETAO,ETALIM,XSI,XSF  AMB1123 

COMMON  /GRP/DMINV,MHINV(101),HMV(101)  AMB1124 

COMMON  /IGRP/KHM  AMB1125 

C  COMPUTE  H(M)  BY  INTERPOLATION  AMB1126 

MEXIT=EM1  AMB1127 

MINV=1.D0/M  AMB1128 

I=KHM-IDINT(MINV/DMINV-l.D-9)-l  AMB1129 

IF(I.GE.l.AND.I.LT.KHM)  GO  TO  1  AMB1130 

PRINT  11,I,KHM,M,MEXIT  AMB1131 

II  FORMATC/1X, 'I,KHM,M,MEXIT=,,2I5,2D14.6/)  AMB1132 
CALL  SOF(f2201')  AMB1133 

1     CONTINUE  AMB1134 

F1=(MINV-MHINV(I+1))/DMINV  AMB1135 

F2=1.D0-F1  AMB1136 

IFCF1 .LT.-1.D-9)  CALL  SOFC2210')  AMB1137 

IFCF2.LT. -l.D-9)  CALL  S0F(f2211')  AMB1138 

H=F1XHMV(I)+F2XHMV(I+1)  AMB1139 

RETURN  AMB1140 

END  AMB1141 

SUBROUTINE  MATCHCT, PSI 0 , MAB, TETAAB)  AMB1142 

C   SUBROUTINE  NUMBER   23  AMB1143 

IMPLICIT  REAL*8(A-H,0-Z,$)  AMB1144 

REALX8  M, MOB, MEXIT, MAB, LAMDOB, KAPAOB  AMB1145 
COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15,AMB1146 

1  G16,G17,G18,G19,G20  AMB1147 

COMMON  /PAR/CO,  ENO,  EMI,  D,  SIGMA,  TLIM,DRO,  EL  0,Q0,  TO,  FACT,  ALOGF,  AMB1148 

1  DPSIO,DTMAX,DETAO,ETALIM,XSI,XSF  AMB1149 

COMMON  /NPAR/NPHI,IPAR,NP,NR,NX,NXS,NS,NSPEC,NS1,NS2,NTAU0,NETA0,  AMB 11 50 

1  NAMB,NCASE,ICASE,IFAN  AMB1151 

COMMON  /GE0M/APF,PAI,PAI2,W,SW,CW, BETA, SBETA,CBETA, PSI 1,SPSI1,     AMB1152 
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1MB      SCRIPT    Al 

1  CPSI1,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,A0,RF,XF,YF,ZF,AMB1153 

2  PHISOF,PHIF,SPHIF,CPHIF,DYMIN,RMIN,XS,DIST,X0,Y0,Z0,  AMB1154 

3  DYO,DEG,PSIN,ST1,CT1,OMEGX,OMEGY,OMEGZ,XSV(21)  AMB1155 
COMMON  /POINT/XP,YP,XCOR,YCOR  AMB1156 
COMMON  /GRP/DMINV,MHINVC101),HMV(101)  AMB1157 
COMMON  /IGRP/KHM  AMB1158 
MEXIT=EM1  AMB1159 
GO  TO  (101,102),IFAN  AMB1160 

101  CONTINUE  AMB1161 
C   FAN  APPROXIMATED  AS  PLANAR  AMB1162 

MAB=DSQRT(1 .D0+G4/DTAN((PSI0-PSIF)/G5)X*2)  AMB1163 

TETAAB  =  PSI0-DASIN(1.D0/MAB)  AMB1164 

GO  TO  100  AMB1165 

102  CONTINUE  AMB1166 
C   COMPUTE  MAB  FROM  THE  INVERSE  PROBLEM  SOLUTION  AMB1167 

C0TAV=1.D0/DTAN(PSI0)  AMB1168 

EVY=YP*DLOG(YP/YCOR)/(YP-YCOR)-l.D0  AMB1169 

PSIN=PSI0  AMB1170 

DO  1  ITER=1,10  AMB1171 

PSI=PSIN  AMB1172 

M=DSQRT(1.D0+G<4/DTANCCPSI-PSIF)/G5)**2)  AMB1173 

M=DMAX1(M,MEXIT)  AMB1174 

CALL  HINTER(M,HM)  AMB1175 

CALL  MFUNC(M,F,ETALF,TETA)  AMB1176 

GOREM=1.DO+G1*M**2  AMB1177 

GOR=M**2-l.D0  AMB1178 

DELT0B=0.5D0*DSQRT(GOR)*(l.D0/(MEXIT*ETALF)  AMB1179 

1      +DSIN(TETA)/M)/DSIN(PSI)+G15*HM/2.D0  AMB1180 

EPSI0B=DELT0B/DSQRT(GOR)-DSIN(TETA)/(M*DSIN(PSI))  AMB1181 

LAMD0B=EPSI0B-DELT0B*GOREM/(GOR*DSQRT(GOR))  AMB1182 

COTN=COTAV+LAMD0B*EVY/DSIN(PSI)**2  AMB1183 

PSIN  =  PAI2-DATAN(C0TN)  AMB1184 

DPSI=PSIN-PSI  AMB1185 

IF(DABS(DPSI).LT.l.D-6)  GO  TO  11  AMB1186 

I  CONTINUE  AMB1187 
PRINT  12,I,ITER,PSI,PSIN,DPSI,M,XP,YP,T  AMB1188 

12    FORMATC/1X, f I , ITER, PSI , PSIN, DPSI , M, XP, YP, T= »//  AMB1189 

1        1X,2I4,7D11.3/)  AMB1190 

CALL  SOF( »2301')  AMB1191 

II  CONTINUE  AMB1192 
C   USING  M0B=M  AS  COMPUTED  FROM  THE  INVERSE  PROBLEM,  FIND  MAB.  AMB1193 

M0B  =  M  AMB1194 

CALL  MFUNCCM,F,ETALF,TETA)  AMB1195 

PSI=TETA+DASIN(1.D0/M)  AMB1196 

CALL  HINTER(M,HM)  AMB1197 

GOREM=l.D0+Gl*M**2  AMB1198 

GOR=M**2-l.D0  AMB1199 

DELT0B=0.5D0*DSQRT(GOR)*(l.D0/(MEXIT*ETALF)  AMB120  0 

1      +DSIN(TETA)/M)/DSIN(PSI)+G15*HM/2.D0  AMB1201 

F0B  =  (G7*GOREM:>x*G2/M  AMB1202 

FAB=F0B*(YP/YCOR)**DELT0B  AMB1203 

CALL  AREAF(FAB,MAB)  AMB1204 

EPSI0B=DELT0B/DSQRT(GOR)-DSIN(TETA)/(M*DSIN(PSI))  AMB1205 

KAPA0B=1.D0  AMB1206 

IF(DABS(PAI2-TETA) .GT.l.D-8)  AMB1207 

1KAPA0B=DTAN(TETA)*EPSI0B  AMB1208 

COSTAB=DC0S(TETA)*(YP/YC0R)*X(-KAPA0B)  AMB1209 

TETAAB=DACOS(COSTAB)  AMB1210 

100   CONTINUE  AMB1211 

RETURN  AMB1212 

END  AMB1213 

SUBROUTINE  AREAF(F,M)  AMB1214 

C   SUBROUTINE  NUMBER   24  AMB1215 

IMPLICIT  REAL*3(A-H,0-Z,$)  AMB1216 

REAL*8  MEXIT,MIN,M,MHINV  AMB1217 
COMMON  /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15,AMB1218 

1             G16,G17,G18,G19,G20  AMB1219 

COMMON  /PAR/CO, ENO, EMI, D, SIGMA, TLIM,DR0, EL 0,Q0, TO, FACT, ALOGF,  AMB1220 

1            DPSIO,DTMAX,DETAO,ETALIM,XSI,XSF  AMB1221 

COMMON  /GRP/DMINV,MHINV(101),HMV(101)  AMB1222 

COMMON  /IGRP/KHM  AMB1223 

C   COMPUTE  MACH  NUMBER  M  FROM  AREA  RATIO  FUNCTION  F  AMB1224 
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AMB       SCRIPT    Al 

C   F=(C2/(G+1))X(1+(G-1)*MXX2))XX((G+1)/(2*(G-1)))/M  AMB1225 

C   INITIAL  GUESS  IS  MIN  AMB1226 

MEXIT=EM1  AMB1227 

E1=CF*MEXIT)X*(1.D0/G2)/G7  AMB1228 

E2=(E1-1.D0)/G1  AMB1229 

E3=DMAX1(E2,MEXITX*2)  AMB1230 

MIN=DSQRT(E3)  AMB1231 

EMN=MIN  AMB1232 

DO  1  1=1,100  AMB1233 

EMO=EMN  AMB1234 

GOREM=l.D0+Gl*EMOXX2  AMB1235 

GOR=EMO**2-l.D0  AMB1236 

F0=CG7XG0REM)**G2/EM0  AMB1237 

DF=FO-F  AMB1238 

C      PRINT  123,I,EM0,EMN,F0,F,DF,G0R,G0REM  AMB1239 

C123   FORMATC1X, ' I , EMO, EMN, FO, F, DF,GOR, GOREM= », 15, 7D12 . 4)  AMB1240 

DFDM=F'3XGOR/(EMOXGOREM)  AMB1241 

DMN  =  DF/DFDM  AMB1242 

EMN=EMO-DMN  AMB1243 

EPSEM=DABS(DMN/EMN)  AMB1244 

IFCEPSEM.LT. l.D-10)  GO  TO  11  AMB1245 

I  CONTINUE  AMB1246 
CALL  SOFC,2401*)  AMB1247 

II  CONTINUE  AMB1248 
M=EMN  AMB1249 
RETURN  AMB1250 
END  AMB1251 
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